function plucked_string_waves_numerical_normal_modes(m)
%% solving the wave equation numerically.
%% used here to demonstrate eigenvalue Sturm-Liouville problems 

%% applied math 105b, spring 2010.
%% Eli Tziperman
clf

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%xs
%% input argument m needs to be a positive integer
if (nargin==0)
  m=1;
  disp('no argument used, assuming input wavenumber m=1;');
  disp('you can specify it instead as plucked_string_waves_numerical_normal_modes(5)');
else
  if m<1 || round(m)~=m
    disp('sorry, argument must be positive integer')
    close all
    return
  end
end


%% domain is [-L,L]
J=200;
L=1;
dx=2*L/J;
x=dx*[0:J]-L;
dt=0.001;
T=100;
N=floor(T/dt);
c2=2;
nplot=10;

if floor(m/2)==m/2; % even
  u_prev=sin(m*pi*x/(2*L));
else                % odd
  u_prev=cos(m*pi*x/(2*L));
end


% boundary conditions:
u_prev(1)=0; u_prev(J+1)=0;
u_now=u_prev; 
u_next=u_prev;
% integrate in time
for n=1:N
  % go over all spatial grid points:
  for j=2:J
    u_xx=(u_now(j+1)-2*u_now(j)+u_now(j-1))/dx^2;
    u_next(j)=dt^2*c2*u_xx+2*u_now(j)-u_prev(j);
  end
  % plot:
  if floor(n/nplot)*nplot==n
    plot(x,u_next)
    axis([-1 1 -1 1])
    title_string=sprintf('time (n)=%d',n);
    h1=title(title_string);
    h2=xlabel('x');
    h3=ylabel('y');
    set([gca,h1,h2,h3],'fontsize',14);
    pause(0.1)
  end
  % prepare for next time step
  u_prev=u_now;
  u_now=u_next;
end

