solve for mathieu duffing type ODE

[Kollányi Tibor<kollanyi at szgtirix dot szgt dot uni-miskolc dot hu>]

Hello,

I would like solve the next mathieu duffing type ODE:
y''+(lambda-2*mu*cos(2*t))y+y^3=0
But there are some problem. Can everybody help me? Thanks.
The code is :

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

int
func (double t, const double y[], double f[],
     void *params)
{
 double mu = *(double *)params;
 double lambda= *((double *)params+1);   
 double epsilon= *((double *)params+2);   
 f[0] = y[1];
 f[1] = -1.0*(lambda-2.0*mu*cos(2.0*t))*y[0]-epsilon*y[0]*y[0]*y[0];
 return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy,
    double dfdt[], void *params)
{
 double mu = *(double *)params;
 double lambda= *((double *)params+1);   
 double epsilon= *((double *)params+2);
 gsl_matrix_view dfdy_mat
   = gsl_matrix_view_array (dfdy, 2, 2);
 gsl_matrix * m = &dfdy_mat.matrix;
 gsl_matrix_set (m, 0, 0, 0.0);
 gsl_matrix_set (m, 0, 1, 1.0);
 gsl_matrix_set (m, 1, 0, 
-1.0*(lambda-2.0*mu*cos(2.0*t))-3.0*epsilon*y[0]*y[0]);
 gsl_matrix_set (m, 1, 1, 0);
 dfdt[0] = 0.0;
 dfdt[1] = -4.0*mu*y[0]*sin(2.0*t);
 return GSL_SUCCESS;
}

int
main (void)
{
 const gsl_odeiv_step_type * T
   = gsl_odeiv_step_rk8pd;

 gsl_odeiv_step * s
   = gsl_odeiv_step_alloc (T, 2);
 gsl_odeiv_control * c
   = gsl_odeiv_control_y_new (1e-6, 0.0);
 gsl_odeiv_evolve * e
   = gsl_odeiv_evolve_alloc (2);

 double mu = 0.2;
 double lambda=1.0;
 double epsilon=0.5;
 double *par;
 par=(double*)malloc(3*sizeof(double));
 par[0]=mu;
 par[1]=lambda;   
 par[2]=epsilon;
 gsl_odeiv_system sys = {func, jac, 2, &par};

 double t = 0.0, t1 = 1.0;
 double h = 1e-6;
 double y[2] = { 0.3, 0.1 };

 while (t < t1)
   {
     int status = gsl_odeiv_evolve_apply (e, c, s,
                                          &sys,
                                          &t, t1,
                                          &h, y);

     if (status != GSL_SUCCESS)
         break;

     printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
   }

 gsl_odeiv_evolve_free(e);
 gsl_odeiv_control_free(c);
 gsl_odeiv_step_free(s);
 return 0;
}
EXECUTING:
.....
----------------------------------------------
1.00000e-06 3.00000e-01 1.00000e-01
6.00000e-06 3.00001e-01 1.00000e-01
3.10000e-05 3.00003e-01 1.00000e-01
1.56000e-04 3.00016e-01 1.00000e-01
7.81000e-04 3.00078e-01 1.00000e-01
7.81000e-04 3.00078e-01 -4.82170e+09
7.81000e-04 3.00078e-01 -8.26019e+09
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan
7.81000e-04 nan nan