# Coupled second order differential equation with NDSolve

I am trying to solve a 2nd order ODE to reproduce a plot. Here are the equations:

dline1[n_, γ_, α_,
Vz_, μ_, η_, Δ_,
r] = -η (D[y[r], r] + 1/r*D[uup[r], r]) + (Vz - μ)*
uup[r] + α*(1/r*udown[r] +
D[udown[r], r]) - Δ*Exp[I*γ]*udown[r];
dline2[n_, γ_, α_,
Vz_, μ_, η_, Δ_,
r] = -α*(D[uup[r], r]) + Δ*Exp[-I*γ]*
uup[r] - η (z'[r] +
1/r*D[udown[r], r] - (n + 1)^2/(4 r^2)*
udown[r]) + (-Vz - μ)*udown[r];


It works okay for Delta=0 and solving analytically:

solin1 = dline1[1, 0, 1, 1, 0, 1, 0, r] /. {uup[r] -> BesselJ[0, z*r],
D[uup[r], r] -> D[BesselJ[0, z*r], r],
y'[r] -> D[BesselJ[0, z*r], {r, 2}], udown[r] -> BesselJ[1, z*r],
D[udown[r], r] -> D[BesselJ[1, z*r], r],
z'[r] -> D[BesselJ[1, z*r], r]};
solin2 = dline2[1, 0, 1, 1, 0, 1, 0, r] /. {uup[r] -> BesselJ[0, z*r],


But it does not work for nonzero Delta and NDSolve, since it deverges to +infinity for one solution and -infinity for other. I tried to make it uncoupled to see what is going on (alpha=Delta=0). The Second equation is Bessel-like, as expected:

soluncoupled2 =
NDSolve[{dline2[1, 0, 0, 1, 0, 1, 0, r] == 0,
z[r] == D[udown[r], r], z[ϵ] == 0.5,
udown[ϵ] == 0},
{udown}, {r, 40}, Method ->
"Automatic"}
] // Flatten;
Plot[Evaluate[{udown[r]} /. soluncoupled2], {r, 0, 40},
PlotRange -> Automatic,
PlotLegends -> {"\!$$\*SubscriptBox[\(u$$, $$↓$$]\)"}]


but the first one is not!

soluncoupled1 =
NDSolve[{dline1[1, 0, 0, 1, 0, 1, 0, r] == 0, y[r] == D[uup[r], r],
y[ϵ] == 0,
uup[ϵ] == 1},
{uup}, {r, ϵ, 40}, MaxSteps -> Infinity,
Method -> {"Automatic"}] // Flatten;

Plot[Evaluate[{uup[r]} /. soluncoupled1], {r, 0, 40},
PlotRange -> Automatic,
PlotLegends -> {"\!$$\*SubscriptBox[\(u$$, $$↑$$]\)"}]


Any ideas? I tried many options for "Methods". Even making u[40]=v[40]=0 does not work. I also changed the boundaries, but still diverges.

• I forgot to mention that trying to solve it analytically for the uncoupled situation in the first one also returns me an error: DSolve::bvnul: For some branches of the general solution, the given boundary conditions lead to an empty solution. Aug 31, 2020 at 19:42
• and solving it all together numerically gives me NDSolve::bvluc: The equations derived from the boundary conditions are numerically ill-conditioned. The boundary conditions may not be sufficient to uniquely define a solution. If a solution is computed, it may match the boundary conditions poorly. and NDSolve::berr: The scaled boundary value residual error of 2.528793742383372*^11 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found. Aug 31, 2020 at 19:44
• With what, exactly, do you need help? Are you seeking symbolic or numerical solutions, and for what parameters? Sep 1, 2020 at 16:10
• @bbgodfrey hey, I am seeking for numerics. I am trying to reproduce a plot that I know it is going to be Bessel like, decreasing in a way that for u[40]=v[40]=0. But my plots are diverging to +inf for u and -inf for v. I don't know what is going on... Sep 1, 2020 at 22:19
• Thanks for your note. I hope you will become a regular contributor to StackExchange. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Oct 21, 2020 at 0:13

soluncoupled1 can be obtained symbolically:

DSolveValue[{uup[r] - Derivative[1][uup][r]/r - uup''[r] == 0,
uup[0] == 1, uup'[0] == 0}, uup[r], r]
(* BesselJ[0, I r] *)


which is equivalent to BesselI[0, r]:

FullSimplify[BesselJ[0, I r] == BesselI[0, r], r > 0]
(* True *)


Because BesselI[0, r] diverges exponentially for large r`, it is not surprising that numerical solutions do the same.

Of course, it may be possible to find some set of parameters for which the coupled set of equations does not diverge. If so, numerical issues would be likely, because avoiding divergences probably would require precise numerical cancellation of some terms.

• Thank you! It makes sense. Oct 20, 2020 at 14:33