function[t,y]=rungekut_sist(a,b,y0,n)
%Efectua n pasos del metodo de Runge-Kutta
%para aproximar la solucion de un sistema de 
%ecuaciones diferenciales 
%(y1',...,yn')=(f1(t,y1,...,yn),...,fn(t,y1,...,yn)) 
%con t en [a,b]
%y con condicion inicial (y1(a),...,yn(a))=y0
%las funciones (f1,...,fn)=fun se guardan en un fichero
%fun.m segun el sistema que queramos resolver
%
m=length(y0);
h=(b-a)/n; t=a:h:b;
y=zeros(n,m);
y(1,:)=y0;

for k=1:n
   k1=fun(t(k),y(k,:));  tk2=t(k)+h/2;
   k2=fun(tk2,y(k,:)+h/2*k1);
   k3=fun(tk2,y(k,:)+h/2*k2);
   k4=fun(t(k+1),y(k,:)+h*k3);
   y(k+1,:)=y(k,:)+h/6*(k1+2*k2+2*k3+k4);
end