%
% Metodo di Runge-Kutta 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]
%
%  Tableau struttura con campi
%   ns     = numero di stadi per il metodo di Runge Kutta
%   a(:,:) = matrice
%   b(:)   = vettore
%   c(:)   = vettore
%
% 
% Uscita
%  x(:)  = vettore soluzione delle 'x'
%  y(:)  = vettore soluzione delle 'y'
%
function [x,y]=RK(f,x0,y0,xf,N,TABLEAU)
  x(1) = x0 ;
  y(1) = y0 ;
  h    = (xf-x0)/N ;
  for k=1:N
    % griglia spazio
    x(k+1) = x(k) + h;
    for s=1:TABLEAU.ns
      % calcolo del coefficente Ks
      ys = y(k) ;
      for j=1:s-1
        ys = ys + TABLEAU.a(s,j)*K(j,:) ;
      end ;
      % xs = x(k) + c(s)*h
      % ys = y(k) + somma( a(s,j)* Kj)
      xs     = x(k)+TABLEAU.c(s)*h ;
      K(s,:) = h*feval(f,xs,ys) ;
    end;
    % y(k+1) = y(k) + somma( b(s) * Ks)
    y(k+1,:) = y(k,:) ;
    for s=1:TABLEAU.ns
      y(k+1,:) = y(k+1,:) + TABLEAU.b(s) * K(s,:) ;
    end ;
  end;
end