function [tout,yout] = rk2(odefun,tspan,y0,options)
%
% Metodo di Eulero modificato
% function [tout,yout] = rk2(odefun,tspan,y0,varargin)
%
% 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;
tout(1) = tspan(1);
n = 0;
yout(:,n+1) = y0(:);
f = zeros(length(y0),2);
while (tout(n+1) < tspan(2))
  n = n+1;
  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));
  yout(:,n+1) = yout(:,n)+k*(b(1)*f(:,1)+b(2)*f(:,2));
  tout(n+1) = tout(n)+k;
end
yout = yout.';