Sample code 6.1.0. Evolve the wave equation based on . The initial profile is

(*Initialisation*)
c=300.0;(*in unit of meter/second*)
Deltax=0.005;(*in unit of m*)
Deltat=Deltax/c;
r=c*Deltat/Deltax;(*Note that r must be less than 1.0 or the wave is numerically unstable*)

xlast=1.0;(* in unit of meter*)
xini=0;(*in unit of meter*)
iini=0;ilast=IntegerPart[(xlast-xini)/Deltax];
tlast=0.01;(* in unit of seconds*)
tini=0;

nini=0;nlast=IntegerPart[(tlast-tini)/Deltat];
x[iini]=xini;x[ilast]=xlast;


(*Initial profile*)
k=1000.0;(*in unit of m^(-2)*)
x0=0.3 (*in unit of meter*);

For[i=iini,i≤ilast,i++,y[i,0]=0.0;y[i,-1]=0;];(*end for i*)

For[i=iini,i≤ilast,i++,
x[i]=x[iini]+i*Deltax;
y[i,0]=Exp[-k(x[i]-x0)^2];
y[i,-1]=Exp[-k((x[i]-Deltax)-x0)^2];
(*y[i,0]=Sin[x[i]/ilast];*)
(*y[i,-1]=0.0;*)
];(*end for i*)
amplitude=Max[Table[Abs[y[i,0]],{i,iini,ilast}]];
txy[0]=Table[{x[idum],y[idum,0]},{idum,iini,ilast}];
ListPlot[txy[0],PlotRange->{{xini,xlast},{-amplitude,amplitude}}]
(*end of initial profile*)


(**Begin time stepping**)
For[n=nini,n≤nlast,n++,

For[i=iini+1,i≤ilast-1,i++,

y[i,n+1]=2(1-(r^2))*y[i,n]-y[i,n-1]+(r^2)*(y[i+1,n]+y[i-1,n]);

(*If[i==IntegerPart[(ilast-iini)/2],
y[i,n+1]=Sin[omega*(n*Deltat)]];(*end if *)*)

(*y[i,n+1]=r^2*(y[i+1,n]-2y[i,n]+y[i-1,n])+2y[i,n]-y[i,n-1];*)
(*Print["i=",i," n=",n,"  ",y[i,n+1]];*)
];(*end for i*)
(*Boundary conditions*)
y[iini,n+1]=0;
y[iini,n+1]=y[iini+1,n];
(*y[ilast,n+1]=y[ilast,n+1];*)
y[ilast,n+1]=y[ilast-1,n];
(*Print[i, " ",n+1," ",y[i,n+1]];*)
(*Boundary conditions*)

txy[n+1]=Table[{x[idum],y[idum,n+1]},{idum,iini,ilast}];
];(*end for n*)
(**End of time stepping**)


(**Visualisation**)
Manipulate[
ListPlot[txy[n],PlotRange->{{xini,xlast},{-2*amplitude,
2amplitude}}]
,{n,nini,nlast-1,1}
](*end Manipulate*)

Sample code 6.1.1. Evolve the wave equation on a composite string, in which speed of the wave travels  c1 on the left half, and c2 on the right half of the string.

(*Initialisation*)
c1=300;(* in unit of m/s*)
c2=150;(* in unit of m/s*)

Deltax=0.005;(*in unit of m*)
Deltat=Deltax/c1;
(*r=c*Deltat/Deltax;*)

xlast=1.0;(* in unit of meter*)
xini=0;(*in unit of meter*)
iini=0;
ilast=IntegerPart[(xlast-xini)/Deltax];
tlast=0.01;(* in unit of seconds*)
tini=0;

nini=0;nlast=IntegerPart[(tlast-tini)/Deltat];
x[iini]=xini;x[ilast]=xlast;


(*Initial profile*)
k=1000.0;(*in unit of m^(-2)*)
x0=0.3 (*in unit of meter*);

For[i=iini,i≤ilast,i++,y[i,0]=0.0;y[i,-1]=0;];(*end for i*)

For[i=iini,i≤ilast,i++,
x[i]=x[iini]+i*Deltax;
y[i,0]=Exp[-k(x[i]-x0)^2];
y[i,-1]=Exp[-k(x[i]-x0)^2];
(*y[i,0]=Sin[x[i]/ilast];*)
(*y[i,-1]=0.0;*)
];(*end for i*)
amplitude=Max[Table[Abs[y[i,0]],{i,iini,ilast}]];
txy[0]=Table[{x[idum],y[idum,0]},{idum,iini,ilast}];
ListPlot[txy[0],PlotRange->{{xini,xlast},{-amplitude,amplitude}}];
(*end of initial profile*)


ihalf=IntegerPart[(ilast-iini)/2];
(**Begin time stepping**)
For[n=nini,n≤nlast,n++,
For[i=iini+1,i≤ilast-1,i++,
If[i≤ihalf,r=c1*Deltat/Deltax;];
If[i>ihalf,r=c2*Deltat/Deltax;];
y[i,n+1]=2(1-(r^2))*y[i,n]-y[i,n-1]+(r^2)*(y[i+1,n]+y[i-1,n]);
(*If[i==IntegerPart[(ilast-iini)/2],
y[i,n+1]=Sin[omega*(n*Deltat)]];(*end if *)*)
];(*end for i*)
(*Boundary conditions*)
y[iini,n+1]=0;y[ilast,n+1]=0;
(*Print[i, " ",n+1," ",y[i,n+1]];*)
(*Boundary conditions*)


txy[n+1]=Table[{x[idum],y[idum,n+1]},{idum,iini,ilast}];
];(*end for n*)
(**End of time stepping**)


(**Visualisation**)
Manipulate[
ListPlot[txy[n],PlotRange->{{xini,xlast},{-amplitude,
amplitude}}]
,{n,nini,nlast-1,1}
](*end Manipulate*)