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.

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

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