function [t,y] = rungekut(fun,t,y0)

% Método de Runge-Kutta (para resolver ecuaciones diferenciales)
% y'(t) = fun(t,y(t)), con y(a)=y0
% t = a, t2, t3, ...,t(n+1)=b
% y: solución en t
m = length(y0);
h = t(2)-t(1);
n = length(t) - 1;
y = zeros(n+1,m);
y(1,:) = y0;

for k = 1:n
    m1 = feval(fun,t(k),y(k,:))';
    m2 = feval(fun,t(k)+h/2,y(k,:)+h/2*m1)';
    m3 = feval(fun,t(k)+h/2,y(k,:)+h/2*m2)';
    m4 = feval(fun,t(k+1),y(k,:)+h*m3)';
y(k+1,:) = y(k,:) + h/6*(m1+2*m2+2*m3+m4);
end