function [t,y]=rungekut(a,b,y0,n,f)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Efectua n pasos del metodo de Runge-Kutta
%para aproximar la solucion de una
%ecuacion diferencial 
%y'=f(t,y) 
%con t en [a,b]y con condicion inicial y(a)=y0
%la funcion f se introduce en modo simbólico 'f' y con variables (t,y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h=(b-a)/n; t=a:h:b;
y=zeros(size(t));
y(1)=y0;

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