function yarr = rungekutta(order, dydt, t0, y0, tfin, deltat)
%% Runge-Kutta 3rd or 4th Order
% Solves first order differential equations with Runge-Kutta methods
% Please predefine t and y as symbolics by typing
% syms t y
% INPUTS:
% ===========
% order (3 is RK3, 4 is RK4)
% dydt symbolic function (differential equation)
% t0 initial time
% y0 initial value
% tfin final time
% deltat time step
% OUTPUTS:
% ===========
% yarr estimated solution at every time step
syms t y

tarr = [t0:deltat:tfin];
yarr = ones(1,length(tarr));
yarr(1) = y0;

if order == 3, % RK3
    for i=1:length(tarr)-1,
        K1 = subs(dydt, [t y], [tarr(i), yarr(i)]);
        K2 = subs(dydt, [t y], [tarr(i)+deltat./2, yarr(i)+K1.*deltat./2]);
        K3 = subs(dydt, [t y], [tarr(i)+deltat, yarr(i)-K1.*deltat+2.*K2.*deltat]);
        yarr(i+1) = yarr(i) + 1./6 .* (K1+4.*K2+K3).*deltat;
    end
else % RK4
    for i=1:length(tarr)-1,
        K1 = subs(dydt, [t y], [tarr(i), yarr(i)]);
        K2 = subs(dydt, [t y], [tarr(i)+deltat./2, yarr(i)+K1.*deltat./2]);
        K3 = subs(dydt, [t y], [tarr(i)+deltat./2, yarr(i)+K2.*deltat./2]);
        K4 = subs(dydt, [t y], [tarr(i)+deltat, yarr(i)+K3.*deltat]);
        yarr(i+1) = yarr(i) + 1./6 .* (K1 + 2.*K2 + 2.*K3 + K4) .* deltat;
    end
end