function [tout,yout] = rk2(odefun,tspan,y0,options)
%
% Metodo di Eulero modificato
% function [tout,yout] = rk2(odefun,tspan,y0,options)
%
% tspan = [t0,t0+tstar]
a(2,1) = 1/2;
c(2) = 1/2;
b(2) = 1;
InitialStep = 0.05;
if (nargin == 4)
  if (isfield(options,'InitialStep'))
    InitialStep = options.InitialStep;
  end
end
k = InitialStep;
n = 1;
tout(n) = tspan(1);
yout(:,n) = y0(:);
xi = zeros(length(y0),2);
%f = zeros(length(y0),2);
while (tspan(2)-tout(n) > eps)
  k = min(k,tspan(2)-tout(n));
  xi(:,1) = yout(:,n);
  xi(:,2) = yout(:,n)+a(2,1)*k*odefun(tout(n)+c(1)*k,xi(:,1));
  yout(:,n+1) = yout(:,n)+k*(b(1)*odefun(tout(n)+c(1)*k,xi(:,1))+...
			     b(2)*odefun(tout(n)+c(2)*k,xi(:,2)));
%  f(:,1) = odefun(tout(n)+c(1)*k,yout(n));
%  f(:,2) = odefun(tout(n)+c(2)*k,yout(n)+a(2,1)*k*f(:,1));
%  yout(:,n+1) = yout(:,n)+k*(b(1)*f(:,1)+b(2)*f(:,2));
  tout(n+1) = tout(n)+k;
  n = n+1;
end
tout = tout.';
yout = yout.';
