%# Number of chunks
N=256;

%# Set up x axis
x=[0:N-1]/N; dx=x(2)-x(1);

%# Set up time step
dt=0.01;

%# Set mass per unit length to one everywhere
mu=ones(1,N);

%# Compute tension which makes wave speed of 1/25 grid point per time step
v=(dx/10)/dt; tau=mu*v^2;

%# Modify mass density from 0.5 to 0.75
mu(0.5*N:0.75*N)=mu(0.5*N:0.75*N)*9;

%# Original function shape, f(x)
y=1.1*exp(-(x-0.25).^2*400).*(x-0.23);
#y=0.04*(1-abs((x-0.25)*8)); y=y.*(y>0);

%# Previous function is just y shifted back by one timestep
yo=1.1*exp(-(x+dt*v-0.25).^2*400).*(x+dt*v-0.23);
#yo=0.04*(1-abs((x+dt*v-0.25)*8)); y=y.*(y>0);

%# Set up graphics
gset yrange [-0.07:0.07];
gset data style linespoints
gset grid

%# Do time steps "forever"
for step=1:1000000
  %# Algorithm to get (d^2 y)/(dx^2)
  yxx=(shift(y,-1)-2*y+shift(y,1))/dx^2;

  %# WAVE EQUATION: (d^2 y)/(dt^2)=(tau/mu) (d^2 y/dx^2)
  ytt=(tau./mu).*yxx;

  %# Algorithm to find next y given current y and (d^2 y)/(dt^2)
  yn=2*y-yo+dt^2*ytt;


  %# Set BOUNDARY CONDITIONS...

  %# Fixed end (zero value)
#  yn(0.5*N)=0; 

  %# Free end (zero slope)
#  yn(0.5*N+1)=yn(0.5*N); 
#  yn(0.5*N+2)=0.; %# Prevents "ghost" wave

  %# Rest of algorithm
  yo=y; y=yn;

  %# Plot every 25th step
  if rem(step,10)==0
    plot(x,y); %# Entire data set
    pause(0);
  endif
end