function [t, x] = RK(fname, tspan, x0, h, G, beta)
% Approssima una soluzione di x'=f(t,x), x(t0)=x0
% nell'intervallo tspan=[tspan(1) tspan(2)].
% con il metodo di RK esplicito a m stadi, definito dalla
% matrice mxm G, e dal vettore beta.
%
% Input:
% fname: nome della funzione f(t,x):RxR^n ----> R^n
% x0: valore iniziale, x0=x(t0)
%  
% h: lunghezza del passo
% G: matrice "Gamma" che definisce il metodo di RK
% beta: vettore "beta" che definisce il metodo di RK
%
% Output:
% t, x: array tali che x(i) e' un'approssimazione di x(t(i)) e
% t(i) = t0 + (i-1)h, i=1:n+1
%
% Esempio:
%
% [t4,x4]=RK(@f3, [0 1], 1/3, 1.e-2, GRK4, bRK4);
% [t3,x3]=RK(@f3, [0 1], 1/3, 1.e-2, GRK3, bRK3);
% [t2,x2]=RK(@f3, [0 1], 1/3, 1.e-2, GHeun, bHeun);
% [t1,x1]=RK(@f3, [0 1], 1/3, 1.e-2, GEulero, bEulero);
%
% x=sol3(t1);
% plot(t1,x,t1,x1,t2,x2,t3,x3,t4,x4)
% legend('esatta','eul','Heun','RK3','RK4')
% Si osserva che RK4 e' quasi coincidente con la sol esatta
% Poi
% err1=abs(x1-x);
% err2=abs(x2-x);
% err3=abs(x3-x);
% err4=abs(x4-x);
% semilogy(t1,err1,t2,err2,t3,err3,t4,err4)
% legend('eul','Heun','RK3','RK4')

   T=tspan(2);t0=tspan(1);
   ro = sum(G, 2);
   npassi = floor( (T-t0)/h );
   tc = t0;
   xc = x0';
   t = tc;
   x = xc;
   fc = feval(fname,tc,xc);
   for j=1:npassi
	  [tc,xc,fc] = RKstep(fname,tc,xc,fc,h,G,ro,beta);
      t = [t;tc];
      x = [x; xc];
   end
