(*constants*)
alpha=1;
x0=0;
x1=1;
t0=0;
t1=2;
dx=.001;
dt=.001;

(*derived constants*)
lambda=alpha*dt/dx;
xlength=IntegerPart[(x1-x0)/dx+1];
tlength=IntegerPart[(t1-t0)/dt+1];

(*helper functions*)
indexToDomain[n_]:=Chop[N[x0+(n-1)*dx],.001];
domainToIndex[x_]:=IntegerPart[(x-x0)/dx+1];
indexToTime[n_]:=Chop[N[t0+(n-1)*dt],.001];
timeToIndex[t_]:=IntegerPart[(t-t0)/dt+1];

(*state vector and past vector*)
u0=ConstantArray[0,xlength];
u1=u0;

(*initial value function and domain*)
f[x_]:=Sin[4*Pi*x]^3;
initialLeft=0;
initialRight=.25;

For[n=domainToIndex[initialLeft],n≤domainToIndex[initialRight],n++,
u0[[n]]=Chop[N[f[indexToDomain[n]]],.001];
]

(*leftC and rightC are the left and right endpoint constants for the boundary derivatives*)
leftC=u0[[1]];
rightC=u0[[xlength]];

(*endpoints*)
b=ConstantArray[0,xlength];
b[[1]]=leftC*lambda^2;
b[[xlength]]=rightC*lambda^2;

(*state change*)
A=SparseArray[{Band[{1,1}]->2(1-lambda^2),Band[{1,2}]->lambda^2,Band[{2,1}]->lambda^2},{xlength,xlength}];

(*initial velocity (i.e. u1)*)
g[x_]:=0;
du=Table[g[indexToDomain[i]],{i,1,xlength}];
u1=.5*A.u0+dt*du+.5*b;

(*solution building*)
answer={u0};
u2=ConstantArray[0,xlength];
For[t=t0,t≤t1,t+=dt,
(*update energy at endpoint*)
leftC=u1[[1]];
rightC=u1[[xlength]];

(*update state vector*)
u2=Chop[A.u1-u0+b,.001];
answer=Join[answer,{u2}];

(*update past vectors*)
u0=u1;
u1=u2;
];

(*how many frames to create*)
frameres=200;
ListAnimate[Table[ListPlot[answer[[i*IntegerPart[(tlength-1)/frameres]+1]],PlotRange->{-1,1}],{i,0,frameres}],40]

Export["1dwave.avi",%]