%% RK method
% Solve the IVP y'(x)=f(x,y) with y(x0)=y0. The function f is given as a
% symbolic function.
% Given: step size h and number of steps. 

function [xRK,yRK] = RK4(f,x0,y0,h,N)

if length(h) == 1
    h = h*ones(1,N);
end
xRK = zeros(1,N+1); yRK = zeros(1,N+1); % Initialization
xRK(1) = x0; yRK(1) = y0;               % Initial condition

for i=1:N                               % For loop over the steps
    k1 = f(xRK(i),yRK(i));              % k1
    k2 = f(xRK(i)+h(i)/2,yRK(i)+h(i)*k1/2);   % k2       
    k3 = f(xRK(i)+h(i)/2,yRK(i)+h(i)*k2/2);   % k3
    k4 = f(xRK(i)+h(i),yRK(i)+h(i)*k3);       % k4
    k = (k1+2*k2+2*k3+k4)/6;            % k
    xRK(i+1) = xRK(i)+h(i);                % Next x
    yRK(i+1) = yRK(i)+h(i)*k;              % Runge Kutta iteration - next y
end