%============ Simulation for a vibrating string ============
clear all;
%---- set dx
dx=10;
dx2=dx*dx;
xsteps=40;
%---- create an x vector whose components are the discrete values of x
x=[0:1:xsteps]*dx;
%---- get the minimum and maximum value of x
minx=min(x);maxx=max(x);
%---- get the size of the x vector
sz0=size(x);
sz=sz0(2);
u=zeros(sz0);
%---- dudt is du/dt
dudt=u;

%---- set the initial condition
%---- you may also change the dt here.
button=questdlg('Choose the initial condition:','Initial condition',' kink ',' waves ',' kink ');
switch button
case ' kink ',
    u(xsteps/2-4)=0.01;u(xsteps/2-3)=0.02;u(xsteps/2-2)=0.03;
    u(xsteps/2-1)=0.04;u(xsteps/2)=0.05;
    u(xsteps/2+1)=0.04;u(xsteps/2+2)=0.03;u(xsteps/2+3)=0.02;
    u(xsteps/2+4)=0.01;
    dt=0.2;
case ' waves '
    u=sin(x*2*pi*5/xsteps/dx)*0.05;
    dt=0.5;
end;

dt2=dt*dt;
tsteps=xsteps*dx/dt;
t=[1:tsteps]*dt;
%---- calculate u(x) for each time steps
for j=1:tsteps,
    %---- plot the function u(x)
    plot(x,u);
    %---- set the ploting range
    axis([minx maxx -0.07 0.07]);
    %---- save the current plot into a movie frame
    m(j)=getframe;
    %---- at each time step, calculate u(x)
    for i=2:sz-1,
        im1=i-1;
        ip1=i+1;
        %---- d2udx2 is (d^2 u)/(dx^2)
        d2udx2=(u(ip1)-2*u(i)+u(im1))/dx2;
%        u(i)=u(i)+(dudt(i)+d2udx2*dt/2)*dt+d2udx2*dt2/2;
        u(i)=u(i)+dudt(i)*dt+d2udx2*dt2/2;
        dudt(i)=dudt(i)+d2udx2*dt;
    end;
%    pause;
end;
%---- save the final movie frame
m(j)=getframe;
%---- play the saved movie
movie(m);
