Find correct boundary value to solve an ODE

This is the first time I am seeking help on this site. I am not sure whether the problem can be sovled. Actually, I am trying to solve an ODE

V[x_] = Q*(1/x0 - 1/x);
g[x_] = 1 - 2*M/x + Q^2/x^2;
U[y_] = y^2 - y^4 + 9/32*y^6;
x0 = 0.15*(1 + 10^-6);
M = 0.102;
Q = 0.09;
e = 0.75;

fEq = f''[x] + (2/x + g'[x]/g[x])*f'[x] + (e^2*V[x]^2)/g[x]^2*f[x] - 1/(2*g[x])*U'[f[x]] // Simplify;

fs[x_] = f0 + x0^3/(2*(x0^2 - Q^2))*U[f0]*(x - x0); (*a function the boundary conditions depends*)

fsol[f0_?NumericQ] := Module[{},
ss = f[x] /.NDSolve[{fEq == 0, f[x0] == f0, f'[x0] == (x0^3 U[f0])/(2 (-Q^2 + x0^2))}, {f[x], f'[x]}, {x, x0, 50}, Method ->"StiffnessSwitching",
Method -> {"ExplicitRungeKutta", Automatic}}][[1]];
ss /. {x -> 50}
](*given f0, NDSolve the OED to get the function value at the other boundary x=50*)

f0sol = f0 /.FindRoot[fsol[f0] == 0, {f0,1.3}] (*search the correct f0 so that f[50]=0*)


the function value at one boundary f[x0]=f0, where f0 is an unkown parameter. I need to search the correct f0 so that the function value at the other end satisfies f[50]=0.

I have tried many times, but failed to get the correct solution. The equation seems very sensitive to the parameter f0, so I wonder if there exists stable solution with correct f0. Is there any method to solve this problem?

• Check the documentation here. By the way, your definition of V[r_] does not contain r. Jul 22 '21 at 3:24
• @bbgodfrey Thanks a lot for pointing out the typo and provding the helpful documentation. Actually, I have tried the shooting method many times, but can not get the right answer. The equation seems very sensitive to f0, so it is hard to search the correct value. I do not know how to solve this problem. Jul 22 '21 at 9:37
• It looks like you're doing something in a Reissner-Nordstrom spacetime. Is x = 50 a special coordinate in your spacetime or is it just a proxy for $r \to \infty$? If it's just a proxy for infinity, can you get the integration to work by picking a smaller value for the outer boundary? Jul 22 '21 at 15:25
• Also, it looks like your inner boundary is just outside the event horizon, and you have used some kind of series extrapolation to approximate the values of f and f' just outside the horizon. What is the actual boundary condition you're trying impose at the horizon? (Otherwise, it seems like f=0 is a perfectly fine solution, but I suspect that would be uninteresting for your purposes.) Jul 22 '21 at 15:46
• @MichaelSeifert Thanks a lot. Yes, I am solving the Klein-Gordon equation on a RN BH background. x=50 is just a proxy for infinity. As you suggested, I will try a smaller value for the outer boundary. The boundary condition on the horizon is just to demand the regularity of the scalar field, so the f and its derivative take finite values there. Jul 22 '21 at 21:42