--- a/scripts/ode/private/runge_kutta_23.m
+++ b/scripts/ode/private/runge_kutta_23.m
@@ -60,10 +60,23 @@
 ## @seealso{runge_kutta_45_dorpri}
 ## @end deftypefn
 
-function [t_next, x_next, x_est, k] = runge_kutta_23 (f, t, x, dt, options = [],
+function [t_next, x_next, x_est, k] = runge_kutta_23 (f, t, x, dt,
+                                                      options = [],
                                                       k_vals = [],
                                                       t_next = t + dt)
 
+  persistent a = [0           0          0;
+                  1/2         0          0;
+                  0           3/4        0];
+  persistent b = [0 1/2 3/4 1];
+  persistent c = [(2/9) (1/3) (4/9)];
+  persistent c_prime = [(7/24) (1/4) (1/3) (1/8)];
+
+  s = t + dt * b;
+  cc = dt * c;
+  aa = dt * a;
+  k = zeros (rows (x), 4);
+
   if (! isempty (options))  # extra arguments for function evaluator
     args = options.funarguments;
   else
@@ -76,34 +89,19 @@
     k(:,1) = feval (f, t, x, args{:});
   endif
 
-  k(:,2) = feval (f, t + (1/2)*dt, x + dt*(1/2)*k(:,1), args{:});
-  k(:,3) = feval (f, t + (3/4)*dt, x + dt*(3/4)*k(:,2), args{:});
-
+  k(:,2) = feval (f, s(2), x + k(:,1) * aa(2, 1).', args{:});
+  k(:,3) = feval (f, s(3), x + k(:,2) * aa(3, 2).', args{:});
+  
   ## compute new time and new values for the unkwnowns
   ## t_next = t + dt;
-  x_next = x + dt.*((2/9)*k(:,1) + (1/3)*k(:,2) + (4/9)*k(:,3)); # 3rd order approximation
+  x_next = x + k(:,1:3) * cc(:); # 3rd order approximation
 
   ## if the estimation of the error is required
   if (nargout >= 3)
     ## new solution to be compared with the previous one
     k(:,4) = feval (f, t_next, x_next, args{:});
-    x_est = x + dt.*((7/24)*k(:,1) + (1/4)*k(:,2) ...
-                     + (1/3)*k(:,3) + (1/8)*k(:,4));
+    cc_prime = dt * c_prime;
+    x_est = x + k * cc_prime(:);
   endif
 
 endfunction
-
-
-%! # We are using the "Van der Pol" implementation.
-%!function [ydot] = fpol (t, y) ## The Van der Pol
-%!  ydot = [y(2); (1 - y(1)^2) * y(2) - y(1)];
-%!endfunction
-%!
-%!shared opts
-%!  opts = odeset;
-%!  opts.funarguments = {};
-%!
-%!test
-%!  [t, y] = runge_kutta_23 (@fpol, 0, [2;0], 0.05, opts);
-%!test
-%!  [t, y, x] = runge_kutta_23 (@fpol, 0, [2;0], 0.1, opts);

--- a/scripts/ode/private/runge_kutta_45_dorpri.m
+++ b/scripts/ode/private/runge_kutta_45_dorpri.m
@@ -67,8 +67,7 @@
                   3/40        9/40       0           0        0          0;
                   44/45      -56/15      32/9        0        0          0;
                   19372/6561 -25360/2187 64448/6561 -212/729  0          0;
-                  9017/3168  -355/33     46732/5247  49/176  -5103/18656 0;
-                  35/384      0          500/1113    125/192 -2187/6784  11/84];
+                  9017/3168  -355/33     46732/5247  49/176  -5103/18656 0];
   persistent b = [0 1/5 3/10 4/5 8/9 1 1];
   persistent c = [(35/384) 0 (500/1113) (125/192) (-2187/6784) (11/84)];
   persistent c_prime = [(5179/57600) 0 (7571/16695) (393/640), ...

--- a/scripts/ode/private/runge_kutta_interpolate.m
+++ b/scripts/ode/private/runge_kutta_interpolate.m
@@ -33,11 +33,10 @@
         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);
+    #{
     case 4
       ## if ode45 is used without local extrapolation this function
       ## doesn't require a new function evaluation.
@@ -85,8 +84,8 @@
 
 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.
+## The function below can be used in an ODE solver to interpolate the
+## solution at the time t_out using 4th order hermite interpolation.
 function x_out = hermite_quartic_interpolation (t, x, der, t_out)
 
   persistent coefs_u_half = ...
@@ -118,3 +117,16 @@
 
 endfunction
 
+## The function below can be used in an ODE solver to interpolate the
+## solution at the time t_out using 3rd order hermite interpolation.
+function x_out = hermite_cubic_interpolation (t, x, der, t_out)
+
+  s = (t_out - t(1)) / (t(2) - t(1));
+  x_out = zeros (size (x, 1), length (t_out));
+
+  for ii = 1:size (x, 1)
+    x_out(ii,:) = (1+2*s).*(1-s).^2*x(ii,1) + s.*(1-s).^2*(t(2)-t(1))*der(ii,1) ...
+                  + (3-2*s).*s.^2*x(ii,2) + (s-1).*s.^2*(t(2)-t(1))*der(ii,2);
+  endfor
+  
+endfunction