I want to solve 2 coupled differential equations with eigenvalue (Ei) like this [Phi]''[r] + (2/r) [Phi]'[r] + mb^2 [Phi][r] + (Ei + g*A[r])^2 [Phi][r] == 0, A''[r] + (2/r) A'[r] + mv^2 A[r] + 2 g (Ei + g*A[r])^2 ([Phi][r])^2 == 0 where mb, mv, and g are constants that we assume 1. The boundary conditions of these equations are [Phi][0] = 1, [Phi]'[0] = 0, A[0] = 0, A'[0] = 0. Because the singularity of r, we assume r = 1e-8. and also I want to find the maximum radius (r) when [Phi][rmax] = 0, [Phi]'[rmax] = 0, A[rmax] = 0, A'[rmax] = 0.

I've tried the code mb = mv = g = 1; b = ParametricNDSolveValue[{[Phi]''[r] + (2/r) [Phi]'[r] + mb^2 [Phi][r] + (Ei + g*A[r])^2 [Phi][r] == 0, A''[r] + (2/r) A'[r] + mv^2 A[r] + 2 g (Ei + g*A[r])^2 ([Phi][r])^2 == 0, [Phi][0.00000001] == 1, [Phi]'[0.00000001] == 0, A[0.00000001] == 0, A'[0.00000001] == 0}, {[Phi], A}, {r, 0.00000001, 100}, {Ei}] but it didn't work when I wanted to find Ei by this val = Map[FindRoot[b[Ei][100], {Ei, #}] &, {1, 2, 3}] and the maximum radius. I want to know where is the problem of my code. I also want to know the plot of Phi and A vs r.