# Solving coupled differential equations with an eigenvalue

I have 2 coupled differential equations with an eigenvalue Ei and want to solve them

 ϕ''[r] + (2/r) ϕ'[r] - mb^2 ϕ[r] + (Ei + g*A[r])^2 ϕ[r] == 0
A''[r] + (2/r) A'[r] - mv^2 A[r] - 2 g (Ei + g*A[r]) (ϕ[r])^2 == 0


where mb, mv and g are constants equal to 1. The boundary conditions of these equations are

ϕ[0] = 1, ϕ'[0] = 0, A[0] = 0, A'[0] = 0


Because of the singularity of r, we assume r = 1*10^-8. What I want to find is the maximum radius (r) when ϕ[rmax] = 0, ϕ'[rmax] = 0, A[rmax] = 0, A'[rmax] = 0. I've with this

mb = mv = g = 1;
b = ParametricNDSolveValue[{ϕ''[r] + (2/r) ϕ'[r] - mb^2 ϕ[r] + (Ei + g A[r])^2 ϕ[r] == 0,
A''[r] + (2/r) A'[r] - mv^2 A[r] - 2 g (Ei + g A[r]) (ϕ[r])^2 == 0,
ϕ[0.00000001] == 1, ϕ'[0.00000001] == 0, A[0.00000001] == 0,
A'[0.00000001] == 0}, {ϕ, A},
{r, 0.00000001, 100}, {Ei}]


but it didn't work when I wanted to find Ei like this

val = Map[FindRoot[b[Ei][100], {Ei, #}] &, {1, 2, 3}]


and the maximum radius. How should I solve this issue ? Also I want to plot ϕ and A vs r.

• What approach do you want to take for finding Ei? b[Ei] actually returns 2 InterpolatingFunction while you only have 1 unknown! If you replace the {ϕ, A} inside ParametricNDSolve with ϕ or A, you'll get a solution, whether the solution is correct or not is another story. BTW, you may want to try this package: library.wolfram.com/infocenter/Demos/4482 Commented Apr 18, 2015 at 4:11
• I want to find Ei where the ϕ and the ϕ' when r (maximum radius) are equal to 0. So I think I can use findroot. Commented Apr 21, 2015 at 7:53
• As mentioned above, ϕ[Ei] == 0 and ϕ'[Ei] == 0 give two equations, use either of them inside FindRoot or RootSearch will give you a result. Commented Apr 21, 2015 at 7:56
• How do I have to write the code of FindRoot like what you've said in Mathematica? Commented Apr 21, 2015 at 9:02
• I noticed that you've modified your equations for several times, are you sure it's correct this time? Maybe you can show us the original problem? Commented Apr 23, 2015 at 3:40

OK, let me explain why your attempt doesn't succeed. Here I'll modify your sample a little for clarity (Of course the troublesome part is leaved untouched ):

lb = 10^-8;
mb = mv = g = 1;

eqn = {ϕ''[r] + (2/r) ϕ'[r] - mb^2 ϕ[r] + (Ei + g*A[r])^2 ϕ[r] == 0,
A''[r] + (2/r) A'[r] - mv^2 A[r] - 2 g (Ei + g*A[r]) (ϕ[r])^2 == 0};
bc = {ϕ[lb] == 1, ϕ'[lb] == 0, A[lb] == 0, A'[lb] == 0};

b = ParametricNDSolveValue[{eqn, bc}, {ϕ, A}, {r, lb, 100}, {Ei}]


As you see, ParametricNDSolveValue returns a ParametricFunction, and again when you give a numeric value to the ParametricFunction, it will return something: it will return what? Let's have a try:

pb = b[1]


Hmm… the calculation stopped at about t = 8.95, far below 100, but it's not the key point here, the key point is, b[1] doesn't return a InterpolatingFunction, it returns a List of InterpolatingFunction! Then what will something like pb[1] return? A expression that FindRoot doesn't understand:

pb[1]


There're many way to circumvent this problem, let me show several here. Since you only have one unknown Ei, only one equation is needed. I'll use ϕ[Ei] for example:

(* Circumvention 1 *)
fb[x_?NumericQ] := b[x][[1]]
FindRoot[fb[Ei][5], {Ei, 1}]

(* Circumvention 2 *)
fb2[x_?NumericQ] := b[x][5] // Through // First
FindRoot[fb2[Ei], {Ei, 1}]

(* Circumvention 3 *)
newb = ParametricNDSolveValue[{eqn, bc}, ϕ, {r, lb, 100}, {Ei}]
FindRoot[newb[Ei][5], {Ei, 1}]

(* Circumvention 4 *)
newb2 = ParametricNDSolve[{eqn, bc}, {ϕ, A}, {r, lb, 100}, {Ei}]
fϕ = ϕ /. newb2
FindRoot[fϕ[Ei][5], {Ei, 1}]


You'll still see some warnings when you execute the above code, but it's mainly because of the nature of your equation and it's another issue, at least FindRoot works this time!

Well, personally I feel the behavior of ParametricNDSolveValue undesirable, why it doesn't return a List of ParametricFunction?

Finally, my intuition told me that you may be interested in this post.

• It's very helpful. Thank you for your help. Commented Apr 28, 2015 at 9:18