2nd order Boundary Value Problem numerical solution (vectors), trying to understand the algorithm

Question

I have a 2nd order BVP:

x''=-Q*x/r 
y''=-Q*y/r 
z''=-Q*z/r 
x(t0)=x0
y(t0)=y0 
z(t0)=z0
x(t1)=x0
y(t1)=y0
z(t1)=z0
r=x^2+y^2+z^2
Q-constant

Need to find x'(t0),y'(t0) and z'(t0)

For solution, I'm choosing the values x'(t0),y'(t0) and z'(t0), and calculate the values for t1. I don't know, how to "calibrate" this values, in case if calculated x(t1),y(t1),z(t1) are wrong.

I know how to implement Runge-Kutta, and read a lot about Shooting method, however I don't know how to change the x',y',z' (3 variables) each time to get the final result.

UPDATE

I use the Shooting method as follows:
y''-P(t)y=F(t) - the form of BVP
y(n+1)=C(n+1)*y(n+2)+D(n+1) 
C(n+1)=1/(2+P(n+1)*h^2-C(n))
D(n+1)=C(n+1)*(D(n)-F(n+1)*h^2)

The problem is, I don't know what are P and F in my case

Answer 1

This is all a lot easier if you generalize a bunch. Given a second order BVP:

$$y''=f(t,y,y') \\ y(a)=y_a \\ y(b)=y_b$$

where $f$ is a given function and $y_a,y_b$ are given vectors. The idea is that you know how to solve the IVP:

$$y''=f(t,y,y') \\ y(a)=y_a \\ y'(a)=y'_a.$$

Therefore, given a value of $y'_a$, you can simulate the IVP (using whatever IVP method you like) and obtain $y(b)$. Let's use $F$ to refer to the mapping $y'_a \mapsto y(b)$; then you just want to solve the equation $F(x)=y_b$ for the vector $x$, which is just a nonlinear system of equations. There are lots of methods out there for nonlinear systems. The only catch here is that the evaluation of $F$ must be treated as a "black box": there is no hope of having an explicit Jacobian of $F$ for general $f$. Still, such methods exist.

This is called the shooting method and it is the simplest method for ODE BVPs. There are faster and more robust methods out there employing different techniques like collocation methods and such.

Answer 2

This is a boundary value problem. No need to supply derivative numbers as initial values. The number of boundary conditions is adequate and correct. Works in all ODE solvers directly. It has a classical closed form solution.

In numerical shoot through, employing sliders for tuning each assumed initial derivative is helpful.

Answer 3

In the ODE system $y''=f(\|y\|)y$ the first derivative $U=\frac{∂y}{∂y'_a}$ satisfies the differential equation $$ U''=f(\|y\|)U+\frac{f'(\|y\|)}{\|y\|^2}yy^TU $$ with $U(a)=0$, $U'(a)=I$. Solving the augmented \begin{align} y''&=f(\|y\|)y,&&y(a)=y_a,\, y'(a)=x\\ U''&=f(\|y\|)U+\frac{f'(\|y\|)}{\|y\|^2}yy^TU,&&U(a)=0,\, U'(a)=I \end{align} system with an ODE -IVP solver will the not only provide the approximation $F(x)=y(b)$, but also the to the same order approximated Jacobian $F'(x)=U(b)$. Then you can apply Newton's method or any more global method based on first order informations.