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

%===============================================================
% Runge-Kutta method for system ordinary differential equations
%===============================================================
clear all
format long

% physical constants/input block
%--------------------------------
% initial conditions
t0 = 0;    % s
x0 = -1;   % m
v0 = 2;    % m/s

% specifications of spring
k  = 1;    % kg/s^2
m  = 2;    % kg

tN = 80;   % end point in s


% define g(x,v,t) and h(x,v,t) specifying the differential equation
%-------------------------------------------------------------------
g  = @(x,v,t)(v);
h  = @(x,v,t)(-k/m * x);


% integrate sstem of ODEs using modified Euler method
%-----------------------------------------------------
[x,v,t] = rk22(g,h,t0,x0,v0,tN,1000);
plot(t,x,'blue'), hold on

[x,v,t] = rk42(g,h,t0,x0,v0,tN,100);
plot(t,x,'red'), hold on

% plot the analytical solution
w   = sqrt(k/m);
A   = sqrt(x0.^2+(v0./w).^2);
phi = atan(w.*x0./v0);
t   = linspace(t0,tN,1000);
plot(t,A*sin(w.*t+phi),'black:')

% label and legend
xlabel('t')
ylabel('x(t)')
title('harmonic oscillator')
legend('1000 2nd order Runge-Kutta steps','100 4th order Runge-Kutta steps','analytical solution','Location','NorthWest')
ylim([-6 6]);


