clear;
format long; format compact

% Recall from Maple
% numbers:={gamma=1, omega=sqrt(5/4), Theta=Pi/4, Omega=2};
gamma=1; omega=sqrt(5/4); Theta=pi/4; Omega=2;
T=15; % Time to stop integrating

% Recall ODE from Maple script:
% diff(theta(t),t,t) + gamma*diff(theta(t),t) + omega^2*sin(theta(t))=0
% Create a function (handle) giving the ODE to be solved:
dy_dt = @(t,y) [y(2); -gamma*y(2) - omega^2 * sin(y(1)) ];

% Now solve the ODE using Matlab's built in function (same as Maple):
[t,Y] = ode45(dy_dt, [0 T], [Theta Omega]);
n_steps = size(t,1)

figure(1); clf
plot(t,Y(:,1),'o--r', t,Y(:,2),'o--b');
title(['Numerical solution']);

figure(2); clf
plot(Y(:,1), Y(:,2), 'o--r');
title(['Phase trajectory']);

%return

% Now let us try to do this ourselves using 100 steps of the Euler method:

N_steps=100; % Number of steps (play around with this)
dt = T/N_steps; % Time step size
t_E = zeros(N_steps+1,1); % Initialize memory for t's
Y_E = zeros(N_steps+1,2); % Initialize memory for Y's
t_E(1)=0; Y_E(1,1)=Theta; Y_E(1,2)=Omega; % Initial conditions
for step=2:N_steps+1
   t_E(step) = t_E(step) + dt; % Time
   Y_E(step,:) = Y_E(step-1,:) + dt * dy_dt(t_E(step-1), Y_E(step-1,:))';
end

figure(1); hold on;
plot(t_E,Y_E(:,1),'sg', t_E,Y_E(:,2),'sk');
legend(['rk45: theta', 'rk45: omega', 'Euler: theta', 'Euler: omega']);

figure(2); hold on;
plot(Y_E(:,1), Y_E(:,2), 's--b');
legend(['rk45', 'Euler']);