rka(x,t,tau,err,derivsRK,param)

function [x, t, tau] = rka(x,t,tau,err,derivsRK,param)
% Adaptive Runge-Kutta routine
% Input arguments -
%   x = current value of dependent variable
%   t = independent variable (usually time)
%   tau = step size (usually timestep) to be attempted
%   err = desired fractional local truncation error
%   derivsrk = right hand side of the ODE; derivsRK is the
%             name of the function with returns dx/dt
%             Calling format derivsRK(x,t,param).
%   param = extra parameters passed to derivsRK
% Output arguments -
%   x = new value of the dependent variable 
%   t = new value of the independent variable
%   tau = suggested step size to be used on next call to rka
tsave = t;  xsave = x;  
safe1 = .9;  safe2 = 4.;  % Safety factors
maxtry = 100;  % Maximum number of attempts
for itry=1:maxtry
  %% Do the two small steps %%
  half_tau = 0.5 * tau;
  xtemp = rk4(xsave,tsave,half_tau,derivsRK,param);
  t = tsave + half_tau;
  x = rk4(xtemp,t,half_tau,derivsRK,param);
  %% Do the single large step %%
  t = tsave + tau;
  xtemp = rk4(xsave,tsave,tau,derivsRK,param);
  %% Compute the error %%
  scale = .5*(abs(x) + abs(xtemp))*err;
  xdiff = x - xtemp;
  errmax = max(abs(xdiff)./(scale + eps));
  %% Estimate new tau value if errmax unacceptable %%
  tau_old = tau;
  tau = safe1*tau_old*errmax^(-0.20);  
  %% Never decrease tau by more than a factor of safe2 %%
  tau = max(tau,tau_old/safe2);
  if errmax < 1,  % If error is acceptable, break out of loop
    break;
  end 
end
%% Never increase tau by more than a factor of safe2 %%
tau = min(tau,safe2*tau_old);