rk4ada2.m

function [t,x,v] = rk4ada2(func, a, b, x0, xd0, h, errmax);
% solution of 2nd order ODE for spring-mass-damper
% using RK 4th order with adaptive step size
%
% USAGE:  [t,x,v] = rk4ada2(func,a,b,x0,xd0,h1,errmax)
%
% input  func    = name of external function to evaluate the RHS
%                  of the ODE for dv/dt (see 'rhs_smd.m')
%        a, b    = limits of integration
%        x0      = initial condition - position
%        xd0     = initial condition - velocity
%        h       = initial stepsize
%        errmax  = desired absolute accuracy
%
% output [t,x,v]  = solution vectors

t = [a];
x = [x0];
v = [xd0];
i = 1;

while t(i) < b

  delta = 2*errmax;
  while abs(delta) > errmax

   % first solution using h1 = h
   h1 = h;
   k1x = v(i);
   k1v = feval(func, t(i)     , x(i)          , v(i)         );
   k2x = v(i) + k1v*h1/2;
   k2v = feval(func, t(i)+h1/2, x(i)+k1x*h1/2 , v(i)+k1v*h1/2);
   k3x = v(i) + k2v*h1/2;
   k3v = feval(func, t(i)+h1/2, x(i)+k2x*h1/2 , v(i)+k2v*h1/2);
   k4x = v(i) + k3v*h1;
   k4v = feval(func, t(i)+h1  , x(i)+k3x*h1   , v(i)+k3v*h1  );
   x1 = x(i) + (k1x + 2*k2x + 2*k3x + k4x)*h1/6;
   v1 = v(i) + (k1v + 2*k2v + 2*k3v + k4v)*h1/6;

   % second solution using h2 = h / 2
   h2 = h / 2;
   k1x = v(i);
   k1v = feval(func, t(i)     , x(i)          , v(i)         );
   k2x = v(i) + k1v*h2/2;
   k2v = feval(func, t(i)+h2/2, x(i)+k1x*h2/2 , v(i)+k1v*h2/2);
   k3x = v(i) + k2v*h2/2;
   k3v = feval(func, t(i)+h2/2, x(i)+k2x*h2/2 , v(i)+k2v*h2/2);
   k4x = v(i) + k3v*h2;
   k4v = feval(func, t(i)+h2  , x(i)+k3x*h2   , v(i)+k3v*h2  );

   % mid-interval values
   tm = t(i) + h2;
   xm = x(i) + (k1x + 2*k2x + 2*k3x + k4x)*h2/6;
   vm = v(i) + (k1v + 2*k2v + 2*k3v + k4v)*h2/6;

   k1x = vm;
   k1v = feval(func, tm     , xm          , vm         );
   k2x = vm + k1v*h2/2;
   k2v = feval(func, tm+h2/2, xm+k1x*h2/2 , vm+k1v*h2/2);
   k3x = vm + k2v*h2/2;
   k3v = feval(func, tm+h2/2, xm+k2x*h2/2 , vm+k2v*h2/2);
   k4x = vm + k3v*h2;
   k4v = feval(func, tm+h2  , xm+k3x*h2   , vm+k3v*h2  );

   x2 = xm + (k1x + 2*k2x + 2*k3x + k4x)*h2/6;
   v2 = vm + (k1v + 2*k2v + 2*k3v + k4v)*h2/6;

   % estimate position error, modify solution, change next step size
   delta = x2 - x1;

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

   if delta ~= 0
      h = h*(errmax/abs(delta))^alpha;
   end

  end  % while delta > errmax

  i = i+1;
  t(i) = t(i-1) + h1;
  x(i) = x2 + delta/15;
  v(i) = v2 + (v2-v1)/15;

end  % while t<b