# Unpacking a Mathematica expression returned by DSolve

I was trying to solve a system of differential equations in Mathematica and had troubles understanding what the solution looked like. So I wanted help to unpack it.

I had a system of two coupled differential equations y1[x] , y2[x] and the solution Mathematica spits out is the following

{{y2 -> Function[{x},
C InverseFunction[
Inactive[Integrate][
1/((Conjugate[K] - C Conjugate[C K]) K^2), {K[
1], 1, #1}] &][-((I x α)/ℏ) + C]],
y1 -> Function[{x},
InverseFunction[
Inactive[Integrate][
1/((Conjugate[K] - C Conjugate[C K]) K^2), {K[
1], 1, #1}] &][-((I x α)/ℏ) + C]]}}


I do not understand function, inverse function, and inactive integrate methods. I was wondering how can I understand the solutions qualitatively and quantitatively.

EDIT : The original code about the two differential equations :

eq1 = I*ℏ*
D[y1[x], x] - α*(Conjugate[y1[x]]*y1[x] - Conjugate[y2[x]]*y2[x])*
y1[x]
eq2 = I*ℏ*
D[y2[x], x] - α*(Conjugate[y1[x]]*y1[x] - Conjugate[y2[x]]*y2[x])*
y2[x]
eqns = {eq1 == 0, eq2 == 0};

• Let us rewrite eq1 = I*\[HBar]*D[y1[x], x] - \[Alpha]*(Abs[y1[x]]^2 - Abs[y2[x]]^2)* y1[x] eq2 = I*\[HBar]*D[y2[x], x] - \[Alpha]*(Abs[y1[x]]^2 - Abs[y2[x]]^2)* y2[x] through Abs, then the solutions are Aug 4, 2021 at 5:06
• {{y2[x] -> C InverseFunction[ Inactive[Integrate][ 1/((Abs[K] - Abs[C K]) (Abs[K] + Abs[C K]) K), {K, 1, #1}] &][-(( I x \[Alpha])/\[HBar]) + C], y1[x] -> InverseFunction[ Inactive[Integrate][ 1/((Abs[K] - Abs[C K]) (Abs[K] + Abs[C K]) K[ 1]), {K, 1, #1}] &][-((I x \[Alpha])/\[HBar]) + C]}} Aug 4, 2021 at 5:11
• I think the system of ODEs under consideration should be solved numerically by ParametricNDSolve. Aug 4, 2021 at 5:16
• In solving DEs, DSolve often has to integrate (obviously) an equation. Sometimes the integral cannot be solved. So it ends up with an equation $F(x,y)=0$ in terms of an integral that is inactivated because it couldn’t be solved. Next DSolve tries to solve for $y$. If the equation is of the form $x = f(y)$, we get $y = f^{-1}(x)$. If $f^{-1}$ cannot be computed (usually can’t when it involves an inactive integral), we $y$ in terms of an inverse function. Aug 4, 2021 at 15:24

In solving DEs, DSolve often has to integrate (obviously) an equation. Sometimes the integral cannot be solved. So it ends up with an equation $$F ( x , y ) = 0$$ in terms of an integral that is inactivated because it couldn’t be solved. Next DSolve tries to solve for $$y$$. If the equation is of the form $$x = f ( y )$$, we get $$y = f^{ − 1} ( x )$$. If $$f^{ − 1}$$ cannot be computed (usually can’t when it involves an inactive integral), we get solution for $$y$$ in terms of an inverse function.

Here's a way to unwind the steps to get the solution in implicit form, still in terms of an unevaluated integral. We'll transform the expressions with three or four replacement rules. The last rule is optional.

(* Convert y -> Function[{x}, a] to y[x] == a *)
{HoldPattern[dv_ -> Verbatim[Function][{iv__}, body_]] :> dv[iv] == body}

(* Solve each equation for the InverseFunction[..][..] in it *)
{eq_Equal :> First@ Solve[eq, Cases[eq, _InverseFunction[_], Infinity]]}

(* Rewrite  $$f^{-1}(a) \rightarrow z$$  as  $$f(z) = a$$  *)
{{Verbatim[InverseFunction][f_][a_] -> z_} :> f[z] == a}

(* Change dummy variable K to something more traditional *)
{K -> ξ}


Apply the rules in succession with ReplaceAll (/.):

sol /.
{HoldPattern[dv_ -> Verbatim[Function][{iv__}, body_]] :>
dv[iv] == body} /.
{eq_Equal :>
First@ Solve[eq, Cases[eq, _InverseFunction[_], Infinity]]} /.
{{Verbatim[InverseFunction][f_][a_] -> z_} :>
f[z] == a} /.
{K -> ξ} Trying to solve a Cauchy problem by

ClearAll["Global*"];
eq1 = I*\[HBar]*D[y1[x],x] - \[Alpha]*(Conjugate[y1[x]]*y1[x] - Conjugate[y2[x]]*y2[x])*y1[x];
eq2 = I*\[HBar]* D[y2[x],x] - \[Alpha]*(Conjugate[y1[x]]*y1[x] - Conjugate[y2[x]]*y2[x])* y2[x];
eqns = {eq1 == 0, eq2 == 0, y1 == I, y2 == 2};
DSolve[eqns, {y1, y2}, x]


, I obtain {} and a warning "DSolve::bvnul: "For some branches of the general solution, the given boundary conditions lead to an empty solution."" The same with other initial values.

The ParametricNDSolve command succeeds:

\[HBar] = 1; \[Alpha] = 2;
eqns = {eq1 == 0, eq2 == 0, y1 == a, y2 == b};
sol = ParametricNDSolve[eqns, {y1, y2}, {x, 1, 2}, {a, b}];
y1 = y1[I, 2] /. sol;
y1[1.5] /. sol


-0.14112 - 0.989992 I`