Solutions of ODE system blow up at certain time step

Question

I am trying to play with the gsl ODE solvers. In particular I am trying to solve a system of two second order ODEs brought into a system of four first order ODEs. I have tried with the implicit Bulirsch-Stoer method. According to GSL website: the method is generally suitable for stiff problems. This stepper requires the Jacobian. I also tried with the explicit embedded Runge-Kutta Prince-Dormand method, which doesn't require the Jacobian.

In the jacobian_hom function I provide both the Jacobian matrix, which takes derivatives wrt. the functions to be solved, and the partial derivatives wrt. the variable t. Bare in mind that t appears both explicitly and within the function Kf, that's why I also give the derivative of Kf, der_Kf.

The system is first solved at initial time ti=1e-5 and final time tf=0.1. Then, I evolve the system in steps of 0.1 for 19 more times, setting ti=tf and tf=tf+0.1. Then I evolve it in steps of 2 for 80 times in the same way. And finally in steps of 10 for 1000 times in the same way. At each step, the system should set the value of the solutions at the previous final time as the new initial conditions (Am I doing it correctly??).

The system should be correct and solvable. My problem is that after some time steps, the solutions blow up and they give -nan after a while, but they should be regular when time becomes large.

I post here an example code. The parameters zM, p, and k that I chose at the beginning of the main can be changed. However, any set of parameters will eventually give me problems with the solutions.

#include <stdio.h>
#include <iostream>
#include <cmath>
#include <cstdlib>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_sf_bessel.h>

using namespace std;

const double gammaE = 0.57721566490153286060651209008240243104215933593992;
const double gW = 0.65407062;
const double gY = 0.35877549;
const double C2 = 3./4.;
const double Y = 1./2.;
const double mW = 1.922745;
const double mB = 0.217089;
const double mphi2 =  0.145381;
const double ml2 = 0.0718046 ;

struct ODEParameters {
    double Meff2;
    double beta;
};

int ode_sys_hom (double t, const double y[], double f[], void *params)
{
  ODEParameters* p = static_cast<ODEParameters*>(params);

  double Meff2 = p -> Meff2 ;
  double beta = p -> beta ;

  double DebW = 1/(2.*M_PI)*( gammaE + log(mW*t/2.) + gsl_sf_bessel_K0(mW*t) );
  double DebB = 1/(2.*M_PI)*( gammaE + log(mB*t/2.) + gsl_sf_bessel_K0(mB*t) );
  double Kf = C2*gW*gW*DebW + Y*Y*gY*gY*DebB; // Kf is a function of t !!

  // System of ODEs
  f[0] = y[2];
  f[1] = y[3];
  f[2] = Meff2 * y[0] - (Kf / beta) * y[1] - (y[2] / t);
  f[3] = Meff2 * y[1] + (Kf / beta) * y[0] - (y[3] / t);

  return GSL_SUCCESS;
}


int jacobian_hom (double t, const double y[], double *dfdy, double dfdt[], void *params)
{

  ODEParameters* p = static_cast<ODEParameters*>(params);
  double Meff2 = p -> Meff2 ;
  double beta = p -> beta ;

  double DebW = 1/(2.*M_PI)*( gammaE + log(mW*t/2.) + gsl_sf_bessel_K0(mW*t) );
  double DebB = 1/(2.*M_PI)*( gammaE + log(mB*t/2.) + gsl_sf_bessel_K0(mB*t) );
  double Kf = C2*gW*gW*DebW + Y*Y*gY*gY*DebB; // Kf is a function of t !!

  double der_DebW = 1/(2.*M_PI)*( 1./t - mW * gsl_sf_bessel_K1(mW*t) );
  double der_DebB = 1/(2.*M_PI)*( 1./t - mB * gsl_sf_bessel_K1(mB*t) );
  double der_Kf = C2*gW*gW*der_DebW + Y*Y*gY*gY*der_DebB; // Kf is a function of t !!

  gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 4, 4);
  gsl_matrix * m = &dfdy_mat.matrix;
  // First row
  gsl_matrix_set (m, 0, 0, 0.0);
  gsl_matrix_set (m, 0, 1, 0.0);
  gsl_matrix_set (m, 0, 2, 1.0);
  gsl_matrix_set (m, 0, 3, 0.0);

  // Second row
  gsl_matrix_set (m, 1, 0, 0.0);
  gsl_matrix_set (m, 1, 1, 0.0);
  gsl_matrix_set (m, 1, 2, 0.0);
  gsl_matrix_set (m, 1, 3, 1.0);

  // Third row
  gsl_matrix_set (m, 2, 0, Meff2);
  gsl_matrix_set (m, 2, 1, - Kf/beta);
  gsl_matrix_set (m, 2, 2, -1./t);
  gsl_matrix_set (m, 2, 3, 0.0);

  // Fourth row
  gsl_matrix_set (m, 3, 0, Kf/beta);
  gsl_matrix_set (m, 3, 1, Meff2);
  gsl_matrix_set (m, 3, 2, 0.0);
  gsl_matrix_set (m, 3, 3, -1./t);

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  dfdt[2] = -( der_Kf/beta ) * y[1] + 1./(t*t) * y[2];
  dfdt[3] =  ( der_Kf/beta ) * y[0] + 1./(t*t) * y[3];

  return GSL_SUCCESS;
}


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

  double zM = 0.1;//stof(argv[1]);
  double p = 1.;//stof(argv[2]);
  double k = 10.;//stof(argv[3]);

  double beta = 0.5 * p / (k-p) / k;
  double alpha = -ml2/2./k + mphi2/2./(k-p) + zM*zM/2./p;
  double Meff2 = alpha / beta;

  cout << "Meff2 = " << Meff2 << "   beta = " << beta << endl;

  // Create an instance of the parameter structure
  ODEParameters params;
  params.Meff2 = Meff2;
  params.beta = beta;

  gsl_odeiv2_system sys_hom_rk = {ode_sys_hom, jacobian_hom, 4, &params}; // Setting the system of ODEs
  gsl_odeiv2_driver * d_hom = gsl_odeiv2_driver_alloc_y_new (&sys_hom_rk, gsl_odeiv2_step_bsimp, 1e-10, 1e-9, 1e-6); // Setting the numerical solver

  double ti = 1e-5;
  double y_hom[4] = { 1.0, 0.0, 0.0, 0.0 };

  int max_steps=0;
  double steps = 0.;

for(int k = 0; k < 3 ; k++){

  if(k==0) {max_steps = 20; steps = 0.1;}
  if(k==1) {max_steps = 80; steps = 2.;}
  if(k==2) {max_steps = 1000; steps = 10.;}

  for (int i = 0; i < max_steps; i++){

    double tf = ti+steps;
      int status_hom = gsl_odeiv2_driver_apply (d_hom, &ti, tf, y_hom);

      if (status_hom != GSL_SUCCESS)
        {
            std::cerr << "Error: GSL integration failed at t = " << ti << std::endl;
            break;
        }

      printf ("Homogenous solutions at t = %.5e, y1 = %.5e, y2 = %.5e, d(y1)/dt = %.5e, d(y2)/dt %.5e\n\n", ti, y_hom[0], y_hom[1], y_hom[2], y_hom[3]);

    ti = tf;
  }
}

  gsl_odeiv2_driver_free (d_hom);

return 0;
}

Any idea on where I am making a mistake? Am I implementing the algorithm correctly?

Thank you very much!