title: Writing Runge-Kutta ODE solver using gsl
body: <p>It's been some time since I did any C/c++, but I wanted to write an ODE solver using the gsl library to solve the following ODE set</p>
  
  <pre><code>$$ u'(r)=up(r)$$
  $$ up'(r)=-(2*(r-1)/(r*(r-2)))*up(r)-((r*r/((r-2)*(r-2)))-(2/r*(r-2)))*u(r) $$
  </code></pre>
  
  <p>so in the gsl notation my y[0]=u, y[1]==up, and the RHS of the above defines f[0] and f[1]. From these definitions one can then compute the Jacobian and dfdr (usually their 'time' variable is called 't' not 'r'). The reason for doing this is because I am having speed issues with Mathematica. I took the gsl example code at the end of their documentation on the ODE solver, and tried to adapt it to my problem as follows:</p>
  
  <pre><code> #include &lt;stdio.h&gt;
   #include &lt;gsl/gsl_errno.h&gt;
   #include &lt;gsl/gsl_matrix.h&gt;
   #include &lt;gsl/gsl_odeiv2.h&gt;
   #include &lt;complex.h&gt;
  
  
   int
   func (double r, const double y[], double f[],
         void *params)
   {
     double mu = *(double *)params;
     f[0] = y[1];
     f[1] = -(2*(r-1)/(r*(r-2)))*y[1]-((r*r/((r-2)*(r-2)))-(2/r*(r-2)))*y[0];
     return GSL_SUCCESS;
   }
  
  /* void tester (double r) {  double outer=-((r*r/((r-2)*(r-2)))-(2/(r*(r-2))));  printf ("%.5e \n", outer); } */      
  
   int
   jac (double r, const double y[], double *dfdy, 
        double dfdt[], void *params)
   {
     double mu = *(double *)params;
     gsl_matrix_view dfdy_mat 
       = gsl_matrix_view_array (dfdy, 2, 2);
     gsl_matrix * m = &amp;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,-((r*r/((r-2)*(r-2)))-(2/(r*(r-2)))));
     gsl_matrix_set (m, 1, 1, -(2*(r-1)/(r*(r-2))));
     dfdt[0] = 0.0;
     dfdt[1] =((1/(r*r))+(1/((r-2)*(r-2))))*y[1]-((4*(1-r)/(r*r*(r-2)*(r-2)))+(4*r/((r-2)*(r-2)*(r-2))))*y[0];
     return GSL_SUCCESS;
   }
  
   int
   main (void)
   {
     /* tester(5);*/
     double om = 2;
     gsl_odeiv2_system sys = {func, jac, 2, &amp;om};
  
     gsl_odeiv2_driver * d = 
       gsl_odeiv2_driver_alloc_y_new (&amp;sys, gsl_odeiv2_step_rk8pd,
                    1e-6, 1e-6, 0.0);
     int i;
     double r = 10, r1 = 100;
     double y[2] = { 0.0000341936, -0.0000572397 };
  
     for (i = 1; i &lt;= 90; i++)
       {
         double ri = 10 + i;
         int status = gsl_odeiv2_driver_apply (d, &amp;r, ri, y);
  
         if (status != GSL_SUCCESS)
      {
        printf ("error, return value=%d\n", status);
        break;
      }
  
         printf ("%.5e %.5e %.5e\n", r, y[0], y[1]);
       }
  
     gsl_odeiv2_driver_free (d);
     return 0;
   }
  </code></pre>
  
  <p>This is giving numbers but they are not the same numbers as Mathematica NDSolve gives, even when I have a low WorkingPrecision and PrecisionGoal. Is there an error in what I've done?</p>