%-----------------------------------------------------------------
% Script 'Discretization Methods', Example 3.7
% (c) 08/2004 by A. Schmidt, IAM, Uni Stuttgart
%-----------------------------------------------------------------
clear all

% Variables
c = 1.0;                            % spring stiffness
m = 1.0;                            % mass
f = 0.0;                            % external force
u0 = 1.0;                           % initial condition displacement
v0 = 0.0;                           % initial condition velocity
t_max = 100.0;                      % time under consideration
dt = 2.001;                         % time step size
n_a = 500;                          % number of sample points 
                                    % (analytical solution)

% Parameters
w = sqrt(c/m);                      % angular frequency

%-----------------------------------------------------------------
% Calculation / time integration by central difference method
t(1) = 0.0;                         % first value (time)
u(1) = u0;                          % initial value (displacement)

u_gp = dt^2/2.0*f/m+(1.0-dt^2/2.0*w^2)*u0-dt*v0;   % ghost point disp. Eq. (3.114)
u(2) = dt^2/m*f+(2.0-dt^2*w^2)*u(1)-u_gp;          % first int. step Eq. (3.113)
t(2) = dt;

for i=2:1:t_max/dt
  t(i+1) = i*dt;
  u(i+1) = dt^2/m*f+(2.0-dt^2*w^2)*u(i)-u(i-1);    % Eq. (3.109)
end

%-----------------------------------------------------------------
% Analytical solution
t_a(1) = 0.0;                       % first value (time)
u_a(1) = u0;                        % initial value (displacement)
for i=1:1:n_a
  t_a(i+1) = i*t_max/n_a;
  u_a(i+1) = u0*cos(w*t_a(i+1));                   % Eq. (3.117)
end

%-----------------------------------------------------------------
% Output / plot
figure(1)
plot(t,u,'r',t_a,u_a,'k')
xlabel('Time / s')
ylabel('Displacement / m')
title('red: central difference method      black: analytical')