Estoy aprendiendo a usar GSL para resolver ODE. Quería resolver el problema del doble péndulo usando las funciones GSL ODE. Decidí usar estas ecuaciones:

(fuente:http://www.physics.usyd.edu.au/~wheat/dpend_html/)

Mi código:

#include <iostream>
#include <cmath>
#include "gsl/gsl_errno.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_odeiv2.h"
#include "constants.h"

double L1;
double L2;
double M1;
double M2;
double T_START;
double T_END;
double S1_ANGLE;
double S2_ANGLE;
double V1_INIT;
double V2_INIT;

int func(double t, const double y[], double f[], void *params) {

    /*
     * y[0] = theta_2
     * y[1] = omega_1
     * y[2] = theta_1
     * y[3] = omega_2
     */

    f[0] = y[1];
    double del = y[2] - y[1];
    double den1 = (M1 + M2) * L1 - M2 * L1 * cos(del) * cos(del);
    f[1] = (M2 * L1 * y[1] * y[1] * sin(del) * cos(del)
            + M2 * G * sin(y[2]) * cos(del) + M2 * L2 * y[3] * y[3] * sin(del)
            - (M1 + M2) * G * sin(y[0])) / den1;

    f[2] = y[3];
    double den2 = (L2 / L1) * den1;
    f[3] = (-M2 * L2 * y[3] * y[3] * sin(del) * cos(del)
            + (M1 + M2) * G * sin(y[0]) * cos(del)
            - (M1 + M2) * L1 * y[1] * y[1] * sin(del)
            - (M1 + M2) * G * sin(y[2])) / den2;

    return GSL_SUCCESS;
}


int main(int argc, char *argv[]) {

    /*
     * Arguments list:
     * 1 - length of pendulum 1
     * 2 - length of pendulum 2
     * 3 - mass of pendulum 1
     * 4 - mass of pendulum 2
     * 5 - start time (seconds)
     * 6 - end time (seconds)
     * 7 - initial angle of 1 pendulum (degrees)
     * 8 - initial angle od 2 pendulum
     * 9 - initial angular velocity of 1 pendulum (deegres per second)
     * 10 - initial angular velocity of 2 pendulum
     */

    if (argc != 11) {
        printf("Wrong number of arguments... \n");
        exit(1);
    }

    //Attribution of arguments
    L1 = atof(argv[1]);
    L2 = atof(argv[2]);
    M1 = atof(argv[3]);
    M2 = atof(argv[4]);
    T_START = atof(argv[5]);
    T_END = atof(argv[6]);
    S1_ANGLE = atof(argv[7]);
    S2_ANGLE = atof(argv[8]);
    V1_INIT = atof(argv[9]);
    V2_INIT = atof(argv[10]);

    //converting to radians
    S1_ANGLE=S1_ANGLE*PI/180.0;
    S2_ANGLE=S2_ANGLE*PI/180.0;
    V1_INIT=V1_INIT*PI/180.0;
    V2_INIT=V2_INIT*PI/180.0;
    printf("L1:%f\nL2: %f\nM1 :%f\nM2:%f\nT_START:%f\nT_END:%f\nS1_ANGLE: %f \nS2_ANGLE: %f\nV1_INIT: %f \nV2_INIT: %f \n",
    L1,L2,M1,M2,T_START,T_END,S1_ANGLE,S2_ANGLE,V1_INIT,V2_INIT);


    gsl_odeiv2_system sys = {func, NULL, 4, NULL};
    gsl_odeiv2_driver *d =
            gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk4, 1e-6, 1e-6, 0.0);



    double y[4] = {S2_ANGLE,V1_INIT,S1_ANGLE,V2_INIT};
    double t = T_START;
    for (int i = 1; i <= 100; i++) {
        double ti = i * (T_END - T_START) / 100.0;
        int status = gsl_odeiv2_driver_apply(d, &t, ti, y);
        printf("%.5e %.5e %.5e %.5e %.5e \n", t, y[0], y[1],y[2],y[3]);
    }


    return 0;
}

Realmente no importa qué parámetros ingrese, siempre obtengo un error:

gsl: driver.c: 354: ERROR: los límites de integración y / o la dirección del paso no son consistentes Se invoca el controlador de error GSL predeterminado.

No entiendo qué causa mi error.