%# 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);

%# Original triangle shape (with a little rounding of peak to help computer)
S=10; u=2*S*(x-0.5); y=(S-(u.*erf(u)+exp(-u.^2)/sqrt(pi)))/S*0.06;

%# Set to zero initial velocity by saying y is the same for previous time step:
yo=y;

%# 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 ends
  yn(1)=0; yn(N)=0;
  %# 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