> restart:with(plots):with(StringTools):

Warning, the name changecoords has been redefined

Warning, the assigned name Group now has a global binding

Set up the parameters: r the interval length l and number of grid point in the spatial direction.

> r:=0.50;l:=1.0;n:=5;

r := .50

l := 1.0

n := 5

Compute mesh spacing h, time step k, and total number of time steps.

> h:=l/n;k:=r*h^2;

h := .2000000000

k := .2000000000e-1

> tsteps:=round(0.3/k);

tsteps := 15

Boundary conditions.

> al:=0.0;ar:=0.0;

al := 0.

ar := 0.

Set up grid xh and initial solution unew from the initial conditions and boundary conditions.

> xh[0+1]:=0.0;

xh[1] := 0.

> unew[0+1]:=al;

unew[1] := 0.

> for i from 1 to n-1 do

> xh[i+1]:=i*h:

> unew[i+1]:=sin(Pi*i*h):

> od:

> xh[n+1]:=l;

xh[6] := 1.0

> unew[n+1]:=ar;

unew[6] := 0.

Graphics: plotting initial condition and initial approximate solution.

> an:=seq([xh[i],unew[i]],i=1..n+1):

> p[1]:=display(pointplot([an],style=line,view=[0..1,0..1],color=blue),
plot(sin(Pi*x),x=0..1),color=red,title=Join(["t=",convert(0,string)],"")):

Time stepping loop.

> for j from 1 to tsteps do

resetting initial solution

> for i from 1 to n+1 do

> uold[i]:=unew[i];

> od;

timestepping the solution (computing solution at next time step)

> unew[0+1]:=al;

> for i from 1 to n-1 do

> unew[i+1]:=evalf(r*uold[i-1+1]+(1-2*r)*uold[i+1]+r*uold[i+1+1]);

> od:

> unew[n+1]:=ar;

plotting exact solution and approximate solution at the new time step

> an:=seq([xh[i],unew[i]],i=1..n+1);

> p[j+1]:=display(pointplot([an],style=line,view=[0..1,0..1],color=blue),plot(exp(-Pi^2*k*j)*sin(Pi*x),x=0..1),color=red,title=Join(["t=",convert(k*j,string)],"")):

> od:

More graphics stuff to create the animation.

> pseq:=seq(p[j],j=1..tsteps+1):

> display(pseq,insequence=true);

[Maple Plot]

>