# HG changeset patch
# User Carlo de Falco <carlo.defalco@polimi.it>
# Date 1449931120 -3600
# Node ID 56d36905893f53bb5469fd5761675af8977c1307
# Parent  f2cd811f0f9ea7914dbd55c047bc13383d9fec7f# Parent  ddc18b909ec77bd7b0376b41198a0965b80d8cf6
codesprint: merge unwanted new head

diff -r f2cd811f0f9e -r 56d36905893f scripts/ode/private/integrate_adaptive.m
--- a/scripts/ode/private/integrate_adaptive.m	Sat Dec 12 09:32:37 2015 -0500
+++ b/scripts/ode/private/integrate_adaptive.m	Sat Dec 12 15:38:40 2015 +0100
@@ -150,8 +150,8 @@
           iout = max (t_caught);
           x(:, t_caught) = ...
             runge_kutta_interpolate (order, [t_old t_new], [x_old x_new],
-                                     tspan(t_caught), new_k_vals, dt,
-                                     options.funarguments{:});
+                                     tspan(t_caught), new_k_vals, dt, func,
+                                     options.funarguments);
 
           istep += 1;
 
diff -r f2cd811f0f9e -r 56d36905893f scripts/ode/private/runge_kutta_interpolate.m
--- a/scripts/ode/private/runge_kutta_interpolate.m	Sat Dec 12 09:32:37 2015 -0500
+++ b/scripts/ode/private/runge_kutta_interpolate.m	Sat Dec 12 15:38:40 2015 +0100
@@ -17,14 +17,15 @@
 ## along with Octave; see the file COPYING.  If not, see
 ## <http://www.gnu.org/licenses/>.
 
-function u_interp = runge_kutta_interpolate (order, z, u, t, k_vals, dt, args)
+function u_interp = runge_kutta_interpolate (order, z, u, t, k_vals, dt, func, args)
 
   switch (order)
 
     ## FIXME: Can interpolations for orders 1-4 simply be deleted? 2015-10-14.
-    #{
+
     case 1
-      u_interp = linear_interpolation (z, u, t);
+      u_interp = interp1 (z, u', t, 'linear');
+
     case 2
       if (! isempty (k_vals))
         der = k_vals(:,1);
@@ -32,6 +33,8 @@
         der = feval (func, z(1) , u(:,1), args);
       endif
       u_interp = quadratic_interpolation (z, u, der, t);
+
+    #{
     case 3
       u_interp = ...
       hermite_cubic_interpolation (z, u, k_vals, t);
@@ -42,7 +45,6 @@
                                        [u(:,i-1) u(:,i)],
                                        k_vals, tspan(counter));
     #}
-
     case 5
       ## ode45 with Dormand-Prince scheme:
       u_interp = hermite_quartic_interpolation (z, u, k_vals, t);
@@ -69,6 +71,19 @@
 
 endfunction
 
+## The function below can be used in an ODE solver to interpolate the solution
+## at the time t_out using 2th order hermite interpolation.
+function x_out = quadratic_interpolation (t, x, der, t_out)
+
+  # coefficients of the parabola
+  a = -(x(:,1) - x(:,2) - der(:).*(t(1)-t(2))) / (t(1) - t(2))^2;
+  b = der(:) - 2*t(1).*a;
+  c = x(:,1) - a*t(1)^2 - b*t(1);
+
+  # evauate in t_out
+  x_out = a*t_out.^2 + b*t_out + c;
+
+endfunction
 
 ## The function below can be used in an ODE solver to interpolate the solution
 ## at the time t_out using 4th order hermite interpolation.
diff -r f2cd811f0f9e -r 56d36905893f scripts/ode/private/starting_stepsize.m
--- a/scripts/ode/private/starting_stepsize.m	Sat Dec 12 09:32:37 2015 -0500
+++ b/scripts/ode/private/starting_stepsize.m	Sat Dec 12 15:38:40 2015 +0100
@@ -43,6 +43,9 @@
 
   ## compute norm of the function evaluated at initial conditions
   y = func (t0, x0);
+  if (iscell (y))
+    y = y{1};
+  endif
   d1 = AbsRel_Norm (y, y, AbsTol, RelTol, normcontrol);
 
   if (d0 < 1e-5 || d1 < 1e-5)
@@ -55,9 +58,12 @@
   x1 = x0 + h0 * y;
 
   ## approximate the derivative norm
+  yh = func (t0+h0, x1);
+  if (iscell (yh))
+    yh = yh{1};
+  endif
   d2 = (1 / h0) * ...
-       AbsRel_Norm (func (t0+h0, x1) - y,
-                    func (t0+h0, x1) - y, AbsTol, RelTol, normcontrol);
+       AbsRel_Norm (yh - y, yh - y, AbsTol, RelTol, normcontrol);
 
   if (max (d1, d2) <= 1e-15)
     h1 = max (1e-6, h0*1e-3);