rk4ada.m

function [t,y] = rk4ada(func, a, b, y0, h1,errmax);
% solution of 1st order ODE using RK 4th order with adaptive step size
%
%      dy/dt = func(t,y)
%
% USAGE:  [t,y] = rk4ada(func,a,b,y0,h1,eps)
%
% input  func    = name of external function to evaluate the RHS
%                  of the ODE
%        a, b    = limits of integration
%        y0      = initial condition
%        h1      = initial stepsize
%        errmax  = desired absolute accuracy
%
% output [t, y]  = solution vectors

t = [a];
y = [y0];
i = 1;

while t(i) < b

   % first solution using h = h1
   h = h1;
   k1 = feval(func, t(i)    , y(i)       );
   k2 = feval(func, t(i)+h/2, y(i)+k1*h/2);
   k3 = feval(func, t(i)+h/2, y(i)+k2*h/2);
   k4 = feval(func, t(i)+h  , y(i)+k3*h  );
   y1 = y(i) + (k1 + 2*k2 + 2*k3 + k4)*h/6;

   % second solution using h = h1 / 2
   h = h1 / 2;
   k1 = feval(func, t(i)    , y(i)       );
   k2 = feval(func, t(i)+h/2, y(i)+k1*h/2);
   k3 = feval(func, t(i)+h/2, y(i)+k2*h/2);
   k4 = feval(func, t(i)+h  , y(i)+k3*h  );
   tm = t(i) + h;
   ym = y(i) + (k1 + 2*k2 + 2*k3 + k4)*h/6;

   k1 = feval(func, tm    , ym       );
   k2 = feval(func, tm+h/2, ym+k1*h/2);
   k3 = feval(func, tm+h/2, ym+k2*h/2);
   k4 = feval(func, tm+h  , ym+k3*h  );
   y2 = ym + (k1 + 2*k2 + 2*k3 + k4)*h/6;

   % estimate error, modify solution, change next step size
   delta = y2 - y1;

   i = i+1;
   t(i) = t(i-1) + h1;
   y(i) = y2 + delta/15;

   if abs(delta) < errmax   % error below tolerance
      alpha = .2;
   else                  % error too big
      alpha = .25;
   end

   h1 = h1*(errmax/abs(delta))^alpha

end