merge changelog

1 files changed, 32 insertions, 30 deletions

diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 5c17836..7f1cf68 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -32,7 +32,7 @@
 #include <gsl/gsl_complex_math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
-#include <gsl/gsl_odeiv.h>
+#include <gsl/gsl_odeiv2.h>
 #include <gsl/gsl_multifit_nlin.h>
 #include <string.h>
 #include <stdio.h>
@@ -1356,38 +1356,38 @@ int ode(int method, double h, double eps_abs, double eps_rel,
         int jac(double, int, const double*, int, int, double*),
         KRVEC(xi), KRVEC(ts), RMAT(sol)) {
 
-    const gsl_odeiv_step_type * T;
+    const gsl_odeiv2_step_type * T;
 
     switch(method) {
-        case 0 : {T = gsl_odeiv_step_rk2; break; }
-        case 1 : {T = gsl_odeiv_step_rk4; break; }
-        case 2 : {T = gsl_odeiv_step_rkf45; break; }
-        case 3 : {T = gsl_odeiv_step_rkck; break; }
-        case 4 : {T = gsl_odeiv_step_rk8pd; break; }
-        case 5 : {T = gsl_odeiv_step_rk2imp; break; }
-        case 6 : {T = gsl_odeiv_step_rk4imp; break; }
-        case 7 : {T = gsl_odeiv_step_bsimp; break; }
-        case 8 : {T = gsl_odeiv_step_gear1; break; }
-        case 9 : {T = gsl_odeiv_step_gear2; break; }
+        case 0 : {T = gsl_odeiv2_step_rk2; break; }
+        case 1 : {T = gsl_odeiv2_step_rk4; break; }
+        case 2 : {T = gsl_odeiv2_step_rkf45; break; }
+        case 3 : {T = gsl_odeiv2_step_rkck; break; }
+        case 4 : {T = gsl_odeiv2_step_rk8pd; break; }
+        case 5 : {T = gsl_odeiv2_step_rk2imp; break; }
+        case 6 : {T = gsl_odeiv2_step_rk4imp; break; }
+        case 7 : {T = gsl_odeiv2_step_bsimp; break; }
+        case 8 : {T = gsl_odeiv2_step_rk1imp; break; }
+        case 9 : {T = gsl_odeiv2_step_msadams; break; }
+        case 10: {T = gsl_odeiv2_step_msbdf; break; }
         default: ERROR(BAD_CODE);
     }
 
-
-    gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin);
-    gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel);
-    gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin);
-
     Tode P;
     P.f = f;
     P.j = jac;
     P.n = xin;
 
-    gsl_odeiv_system sys = {odefunc, odejac, xin, &P};
+    gsl_odeiv2_system sys = {odefunc, odejac, xin, &P};
+
+    gsl_odeiv2_driver * d =
+         gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel);
 
     double t = tsp[0];
 
     double* y = (double*)calloc(xin,sizeof(double));
     int i,j;
+    int status;
     for(i=0; i< xin; i++) {
         y[i] = xip[i];
         solp[i] = xip[i];
@@ -1396,22 +1396,24 @@ int ode(int method, double h, double eps_abs, double eps_rel,
        for (i = 1; i < tsn ; i++)
          {
            double ti = tsp[i];
-           while (t < ti)
-             {
-               gsl_odeiv_evolve_apply (e, c, s,
-                                       &sys,
-                                       &t, ti, &h,
-                                       y);
-               // if (h < hmin) h = hmin;
-             }
+           
+           status = gsl_odeiv2_driver_apply (d, &t, ti, y);
+     
+           if (status != GSL_SUCCESS) {
+                  printf ("error in ode, return value=%d\n", status);
+                  break;
+                }
+
+//           printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
+           
            for(j=0; j<xin; j++) {
                solp[i*xin + j] = y[j];
            }
          }
 
     free(y);
-    gsl_odeiv_evolve_free (e);
-    gsl_odeiv_control_free (c);
-    gsl_odeiv_step_free (s);
-    return 0;
+    gsl_odeiv2_driver_free (d);
+    
+    return status;
 }
+