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)) ];

% 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)

% 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

step=1; % Initialize
t_E(step)=0; Y_E(step,1)=Theta; Y_E(step,2)=Omega; % Initial conditions

% And now do the Euler algorithm usign a for loop
for step=2:N_steps+1
   t_E(step) = t_E(step-1) + dt; % Time
   Y_E(step,:) = Y_E(step-1,:) + dt * dy_dt(t_E(step-1), Y_E(step-1,:))';
end

% Plot the solution
%-----------------------------------
 
if(0)

   figure(1); clf
   plot(t,Y(:,1),'o--r'); hold on
   plot(t,Y(:,2),'s--b');
   title('Numerical solution');
   legend(['theta', 'd(theta)/dt']);

   figure(2); clf
   plot(Y(:,1), Y(:,2), 'o--r');
   xlabel('theta'); ylabel('d(theta)/dt');
   title('Phase trajectory');

else

   figure(1); clf;
   plot(t, Y(:,1), 'o--r'); hold on;
   plot(t_E, Y_E(:,1), 's--b');
   legend(['rk45','Euler']);

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

end
