clear
clf

f=0;

L = 10;
n = 100;
Nstep = 200;
x = linspace(-L,L,n);
velocity = 1;
k = (2*pi/L)/2; %wave number (2pi/wavelength L
p = 1; %period/octave

dx = L/(n-1);
dt = 1/velocity*dx;

A = 10;
%u = A*sin(p*k*x); %sine wave regular standing wave
%u = A*exp(-L*x.^2/L);
%u = A*exp(-x.^4/L);
%u = A*sin(x).*exp(-x.^2/L);
%u = A*cos(x).*exp(-x.^2/L);
%u = A*sin(x).^2.*exp(-x.^2/L);
u = A*cos(x).^2.*exp(-x.^2/L);



title_str= 'Wave equation: u_{tt} -c^2 u_{xx}=0';

v = u; %u(i)
U = u; %u(i-1)
v_initial = u; %constant, original plot

uxx = zeros(1,n); %u(i+1)

for time = 1:Nstep
    t = time*dt;
    y1 = (A/2)*sin((p*k*x-(2*p*k*t)));
    y2 = (A/2)*sin((p*k*x+(2*p*k*t)));
    
    uxx(1) = (v(2)-2*v(1)+v(n))/(dx*dx);
    uxx(2) = (v(3)-2*v(2)+v(n))/(dx*dx);
    uxx(:,n) = (v(1)-2*v(n)+v(n-1))/(dx*dx);
    uxx(:,n-1) = (v(n)-2*v(n-1)+v(n-2))/(dx*dx);
    
    for i=3:n-2
        uxx(:,i) = (v(i+1)-2*v(i)+v(i-1))/(dx*dx);
    end
    
    u = (dt*dt*velocity^2*uxx+2*v)-U;
    
    %Reflecting boundary
    u(1) = 0;
    u(n) = 0;
    
    %update arrays
    U = v; 
    v = u;
    
    %plot(x,y1,'b',x,y2,'g',x,u,'k',-L:L/p:L,0,'or','LineWidth',2);
    plot(x,v_initial,'-.b',x,u,'k','linewidth',3);
    title(title_str)
   
    axis([-L L -1.5*A 1.5*A]);
    drawnow;
    
    
    f=f+1;
    saveas(gcf, ['standing_wave_cos2exp',sprintf('%03d',f), '.png'])
    
end

