function [tout,yout] = rk2(odefun,tspan,y0)
%
% function [tout,yout] = rk2(odefun,tspan,y0)
%
% RK2 Eulero modificato
a(2,1) = 1/2;
c(2) = 1/2;
b(2) = 1;
tout = tspan;
h = diff(tspan);
nmax = length(h);
n = 0;
yout = zeros(size(y0,1),nmax+1);
yout(:,n+1) = y0;
while (n < nmax)
  n = n+1;
  k{1} = odefun(tspan(n),yout(:,n));
  k{2} = odefun(tspan(n)+c(2)*h(n),yout(:,n)+a(2,1)*h(n)*k{1});
  yout(:,n+1) = yout(:,n)+h(n)*(b(1)*k{1}+b(2)*k{2});
end

