function [tout,yout] = rk23(odefun,tspan,y0,options)
%
% function [tout,yout] = rk23(odefun,tspan,y0,options)
%
a(2,1) = 1/2;
a(3,1) = -1;
a(3,2) = 2;
c(2) = 1/2;
c(3) = 1;
b(2) = 1;
bhat(1) = 1/6;
bhat(2) = 2/3;
bhat(3) = 1/6;
InitialStep = 0.01;
RelTol = 1e-3;
AbsTol = 1e-6;
if (nargin == 4)
  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
k = InitialStep;
n = 1;
tout(n) = tspan(1);
yout(:,n) = y0(:);
f = zeros(length(y0),3);
while (tspan(2)-tout(n) > eps)
  k = min(k,tspan(2)-tout(n));
  f(:,1) = odefun(tout(n),yout(:,n));
  f(:,2) = odefun(tout(n)+c(2)*k,yout(:,n)+a(2,1)*k*f(:,1));
% rk2 (Heun)
  yout(:,n+1) = yout(:,n)+k*(b(1)*f(:,1)+b(2)*f(:,2));
  f(:,3) = odefun(tout(n)+c(3)*k,yout(:,n)+a(3,1)*k*f(:,1)+a(3,2)*k*f(:,2));
% prima si calcola la differenza
  err = k*((bhat(1)-b(1))*f(:,1)+(bhat(2)-b(2))*f(:,2)+bhat(3)*f(:,3));
  normyout = norm(yout(:,n+1),Inf);
  normerr = norm(err,Inf);
  if (normerr > AbsTol+normyout*RelTol)
    warning('time step rejected')
  else
% rk3
    yout(:,n+1) = yout(:,n+1)+err;
    tout(n+1) = tout(n)+k;
    n = n+1;
  end
  k = min(2,max(0.6,abs(0.7*...
			((AbsTol+normyout*RelTol)/normerr)^(1/3))))*k;
end
yout = yout.';
