%% This is not quite the wave equation, as I've added some adhoc damping

figure(1); clf;

x = linspace(0, pi, 300);

a = 1;

n = 1:50;   % how many terms in truncated Fourier Series
L = pi/10;  % where to pluck
damping = 0.01;


% copy-paste from CAS: see course on Fourier Series
a1 = (-L .* cos (L .* n) ./ n + sin (L .* n) ./ n .^ 2) ./ L;
a2 = L .* n .* cos (L .* n) ./ (L .* n .^ 2 - pi .* n .^ 2) - pi .* n .* cos (L .* n) ./ (L .* n .^ 2 - pi .* n .^ 2) - sin (L .* n) ./ (L .* n .^ 2 - pi .* n .^ 2);
cn = a1 + a2;

H = plot(x, 0*x, 'linewidth', 3);
ylim([-1.4 1.7])
xlim([0 pi])
xlabel('x');
ylabel('u');
grid on

for i=0:500
  t = 0.1*i;
  %u = an .* sin(a*n*t).*sin(n*x);
  u = zeros(size(x));
  for j=1:length(n)
    % exp(-n^2): an ad-hoc damping term
    u = u + cn(j) .* exp(-damping*n(j)^2*t) .* cos(a*n(j)*t) .* sin(n(j)*x);
  end
  %clf;
  %plot(x, u, 'linewidth', 3);
  set(H, 'ydata', u);
  pause(0.01)
  if (i==0)
    pause
  end
end
