function [tout,yout] = rk23(odefun,tspan,y0,varargin)
%
% function [tout,yout] = rk23(odefun,tspan,y0,varargin)
%
% RK23
a(2,1) = 1/2;
a(3,1) = -1;
a(3,2) = 2;
c(2) = 1/2;
c(3) = 1;
b2(1) = 0;
b2(2) = 1;
b3(1) = 1/6;
b3(2) = 2/3;
b3(3) = 1/6;
n = 0;
yout(:,n+1) = y0;
InitialStep = 0.01;
RelTol = 1e-3;
AbsTol = 1e-6;
if (nargin == 4)
  options = varargin{1};
  if (isfield(options,'RelTol'))
    RelTol = options.RelTol;
  end
  if (isfield(options,'AbsTol'))
    AbsTol = options.AbsTol;
  end
  if (isfield(options,'InitialStep'))
    InitialStep = options.InitialStep;
  end
end
h = InitialStep;
tout(1) = tspan(1);
while (tout(n+1) < tspan(2))
  n = n+1;
  h = min(h,tspan(2)-tout(n));
  k{1} = odefun(tout(n),yout(:,n));
  k{2} = odefun(tout(n)+c(2)*h,yout(:,n)+a(2,1)*h*k{1});
% rk2
  yout(:,n+1) = yout(:,n)+h*(b2(1)*k{1}+b2(2)*k{2});
  k{3} = odefun(tout(n)+c(3)*h,yout(:,n)+a(3,1)*h*k{1}+a(3,2)*h*k{2});
  err = h*((b3(1)-b2(1))*k{1}+(b3(2)-b2(2))*k{2}+b3(3)*k{3});
% rk3
  yout(:,n+1) = yout(:,n+1)+err;
  tout(n+1) = tout(n)+h;
  normyout = norm(yout(:,n+1));
  normerr = norm(err);
  if (normerr > AbsTol+normyout*RelTol)
    disp('time step rejected')
    n = n-1;
  end
  h = min(2,max(0.6,abs(0.9*((AbsTol+normyout*RelTol)/...
			     normerr)^(1/3))))*h;
end
