function [tt,u]=rk2(odefun,tspan,y0,Nh,varargin)
% RK2 Risolve un sistema di e.d.o. con Runge-Kutta2
% cioè a 2 stadi(noto anche come metodo di Heun)
% e passo costante
% b_i=0.5   c_1=0  c_2=1
% INPUT :  nome del file contenente la e.d.o.
%          vettore contenente t0 e tfinale
%          condizione iniziale
%          numero di intervalli
% OUTPUT:  vettore degli istanti temporali
%          vettore con la soluzione approssimata

if nargin==4
    varargin{1}=[] ;
    varargin{2}=[] ;
end

h=(tspan(2)-tspan(1))/Nh; 
tt=tspan(1):h:tspan(2) ;
%tt=linspace(tspan(1),tspan(2),Nh+1)';

y0=y0(:);    % genera sempre un vettore colonna
d=length(y0); % dimensione del sistema   
u=zeros(Nh+1,d);
% inizializzazione del vettore (matrice se d>1) u
u(1,:)=y0' ;
hh=0.5*h ;  %h*b_i
for n=1:Nh
    % u_n+1=u_n+0.5*h*(K1+K2)
    wn=u(n,:) ;
    K1=feval(odefun,tt(n),wn,varargin{1},varargin{2}) ;
    y=wn+h*K1 ;
    K2=feval(odefun,tt(n+1),y,varargin{1},varargin{2}) ;
    w=wn+hh*(K1+K2) ;
    u(n+1,:)=w' ;      
end

plot(tt,u)

if strcmp(odefun,'dinamica_pop')
    title('Dinamica di una popolazione')
end