% Damped Simple Harmonic Oscillator solved by Euler's Method
%
% This routine implements a numerical (Euler's Method) solution to the
% damped simple harmonic oscillator.


% Physical Parameters of the Oscillator
w0=10; beta=1;

w=sqrt(w0^2-(beta/2)^2);
phi=atan(-1/20);
x0=1/cos(phi);

% Parameters for the Integration
h = 0.003;                  % Time step in seconds
duration=15;                % Duration of solution (seconds)
Nstep=round(duration/h);    % Resulting number of steps needed.
t=0:h:Nstep*h;              % Resulting time vector (seconds)

% Reserve some memory for the solution vector (these lines are optional but
% do speed up the code.)
x=zeros(size(t));   % x is position (meters)
y=zeros(size(t));   % y is velocity (meters)

% Set the initial conditions to be the first elements of x and y.
x(1)=x0;             % Initial position (meters)
y(1)=0;             % Initial velocity (meters/second)

% Iterate to find solution with Euler's method
for n=2:Nstep
    
    f=y(n-1);                   % f(x,y,t) is dx/dt
    x(n)=x(n-1)+f*h;            % Use the last position to find the new position

    g=-beta*y(n-1)-w0^2*x(n-1); % g(x,y,t) is dy/dt.
    y(n)=y(n-1)+g*h;            % Use the last velocity to find the new velocity
    
end

% The exact solution (solved on paper) is
xexact=x(1)*exp(-beta/2*t).*cos(w*t+phi);
yexact=x(1)*(-beta/2*exp(-beta/2*t).*cos(w*t+phi)-w*exp(-beta/2*t).*sin(w*t+phi));
GDE=max(xexact-x);

% Compare the exact solution and the numerical solution
figure(1);
plot(t,x,'bo',t,xexact,'r','linewidth',1,'markersize',2); hold on;
title(['GDE = ',num2str(GDE)]);
grid on; hold off;