% Approximates solutions to mx'' + cx' + kx = cos(2 pi fd t) using an
%  explicit time stepping finite difference scheme.  To use a different
%  driving force, pass the vector fj = f(tj) in rather than assigning 
%  fj = cos(2 pi fd t).

clear all
m = 0.01;     % system mass  (kg)
k = 100;      % spring constant  (N/m)
c = 0.1;      % damping coefficient  (N*s/m)
a = 0;        % initial position
b = 0;        % initial velocity
h = 1/44100;  % time step length (s)
Tf = 0.25;    % stop time (s)
fd = 80;      % driving frequency (Hz)

tj = [0:h:Tf];                  % sample times
fj = cos(2*pi*fd*tj);
c1 = 1/((m/h^2) + (c/(2*h)));   % constants for iteration
c2 = (2*m)/h^2 - k;
c3 = -m/h^2 + c/(2*h);

xj(1) = a;                       % set initial conditions x(0) = a,
xj(2) = a + b*h;                 %    x'(0) = b.

for j = 2:length(tj)-1,          % main iteration loop
  xj(j+1) = c1*(fj(j) + c2*xj(j) + c3*xj(j-1));
end

plot(tj,xj)