This is the radial equation for a free particle. To solve:

In[3]:= DSolve[ \[Psi]''[r] + 2/r \[Psi]'[r] + 
   m^2 \[Psi][r] == -E^2 \[Psi][r], \[Psi][r], r]

Out[3]= {{\[Psi][r] -> (E^(-Sqrt[-E^2 - m^2] r) C[1])/r + (
    E^(Sqrt[-E^2 - m^2] r) C[2])/(2 Sqrt[-E^2 - m^2] r)}}

C[2]=0 to satisfy the boundary condition at infinity, C[1] can be determined by normalization.

Prof, now 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][r_max] = 0, \[Phi]'[r_max] = 0, A[r_max] = 0, A'[r_max] = 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.