% Unit 5 (MATLAB Course @ UAM)
% (c) Alexander Knebe

%=====================================================================
% Runge-Kutta method for 2 coupled ordinary differential equations
%=====================================================================
function [x,v,t] = rk42(g, h, t0, x0, v0, tN, N)


% initialisation block
%----------------------
% step size
dt   = (tN-t0)/(N-1);

% arrays
t    = linspace(t0,tN,N);
x    = zeros(1,N);
v    = zeros(1,N);
%t(1) = t0; this is automatically initialized by linspace()
x(1) = x0;
v(1) = v0;

% integration block
%-------------------
for i=1:N-1
    
  
    % trial steps
    k1 = g(x(i),           v(i),          t(i));
    m1 = h(x(i),           v(i),          t(i));
    
    k2 = g(x(i)+k1*dt/2,   v(i)+m1*dt/2,  t(i)+dt/2);
    m2 = h(x(i)+k1*dt/2,   v(i)+m1*dt/2,  t(i)+dt/2);
    
    k3 = g(x(i)+k2*dt/2,   v(i)+m2*dt/2,  t(i)+dt/2);
    m3 = h(x(i)+k2*dt/2,   v(i)+m2*dt/2,  t(i)+dt/2);
    
    k4 = g(x(i)+k3*dt,     v(i)+m3*dt,    t(i+1));
    m4 = h(x(i)+k3*dt,     v(i)+m3*dt,    t(i+1));

        
    % actual steps
    x(i+1) = x(i) + 1/6 * (k1+2*k2+2*k3+k4) * dt;
    v(i+1) = v(i) + 1/6 * (m1+2*m2+2*m3+m4) * dt;


end

