%
% Metodo di Collatz o Eulero Migliorato per sistemi di equazioni differenziali
% 
% Ingresso
%  f      = funzione per il sistema y'=f(x,y)
%  x0, y0 = dato iniziale del problema y(x0) = y0
%  xf     = estremo superiore dell'intervallo di calcolo [x0,xf]
%  N      = numero di intervalli in cui suddividere  [x0,xf]
% 
% Uscita
%  x(:)  = vettore soluzione delle 'x'
%  y(:)  = vettore soluzione delle 'y'
%
function [x,y]=Collatz(f,x0,y0,xf,N)
  x(1)   = x0;
  y(1,:) = y0;
  h      = (xf-x0)/N ;
  for k=1:N
    % griglia spazio
    x(k+1) = x(k) + h;
    % mezzo passo di Eulero
    % y(k+1/2) = y(k) + (h/2)*f(x(k),y(k))
	fk  = feval(f,x(k),y(k,:));
    ykh = y(k,:) + (h/2).*fk;
    % correzione di Collatz
    % y(k+1) = y(k) + h*f(x(k+1/2),y(k+1/2))
	fkh = feval(f,x(k) + h/2,ykh);
    y(k+1,:) = y(k,:) + h.*fkh ;
  end;
end