[Help-gsl] Boundary value problem for system of ODEs

[Ruben Farinelli]

Hi,
I have been working for a long time on a complicate physical problem.
I have a set of ODEs, three of which are first-order, and the fourth
is of second order.
In the latter case unfortunately, I have conditons for its derivative
R'[0]=0,
while the second is actually an asymptotic boundary condition, namely
the function must tend to zero for large value of the variable.
Of course with a change of variable the system becomes N+1
first order ODEs.

Actually I have implemented some kind of shooting method, the main
issue is that the function has an exponential-decay behavior and
the result seems to be rather dependent on the adopted integrator.

Finally I decided to compute the complicate Jacobian to test
the gsl_odeiv_bsimp which I read should be the most powerful.

A-part from any welcome suggestion in approaching BVPs, I have
a doubt. Namely, the right-hand side of the system contains not only
the functions y_i but also their first derivatives, labeled F_i, or y'_i.
I mean something like

F[0]=f_0{y[0]...y[n], F[1].....F[n]}
F[1]=f_1{y[0]...y[n], F[0].....F[n]}
etc etc

The GSL jacobian function arguments are
(double t, const double y[], double *dfdy, double dtdy[], void *params)

but I don't see where the functions derivative y'_i are stored.
They are still present when computing the Jacobian, but apparently
they are not read.

Is there something wrong ?

Thank you for your help !
Ruben