{VERSION 5 0 "IBM INTEL NT" "5.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Times" 1 12 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 
0 0 1 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" 0 21 "" 0 1 0 0 0 1 0 0 0 0 2 0 0 0 0 1 }{CSTYLE "MYTEXT" 
-1 256 "Times" 1 14 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{PSTYLE "Normal" -1 
0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 
0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 1" -1 3 1 {CSTYLE "" -1 -1 "Ti
mes" 1 18 0 128 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 6 6 1 0 1 0 2 2 0 1 }
{PSTYLE "R3 Font 0" -1 256 1 {CSTYLE "" -1 -1 "Courier" 1 12 255 0 0 
1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "R3 Font 2
" -1 257 1 {CSTYLE "" -1 -1 "Courier" 1 12 0 0 255 1 2 2 2 2 2 2 1 1 
1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Normal" -1 258 1 {CSTYLE "
" -1 -1 "Times" 1 12 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }3 1 0 0 0 0 1 0 1 0 
2 2 0 1 }}
{SECT 0 {SECT 0 {PARA 3 "" 0 "" {TEXT -1 37 "Diffusion equation---nume
rical solver" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 62 "This program numer
ically solves a diffusion equation on (0,1):" }}{PARA 258 "" 0 "" 
{XPPEDIT 18 0 "diff(u(x,t),t) = D*diff(u(x,t),`$`(x,2));" "6#/-%%diffG
6$-%\"uG6$%\"xG%\"tGF+*&%\"DG\"\"\"-F%6$-F(6$F*F+-%\"$G6$F*\"\"#F." }
{TEXT -1 1 "+" }{XPPEDIT 18 0 "lambda*f(u(x,t));" "6#*&%'lambdaG\"\"\"
-%\"fG6#-%\"uG6$%\"xG%\"tGF%" }}{PARA 0 "" 0 "" {TEXT -1 70 "with Diri
chlet boundary problem u(0,t)=u(L,t)=0 and initial condition " }
{XPPEDIT 18 0 "u(x,0) = sin(Pi*x);" "6#/-%\"uG6$%\"xG\"\"!-%$sinG6#*&%
#PiG\"\"\"F'F." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}}
{SECT 0 {PARA 3 "" 0 "" {TEXT -1 16 "Matrix Algorithm" }}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 91 "L is the length of spatial interval; k is the time step; \+
h is the grid size(spatial step); " }}{PARA 0 "" 0 "" {TEXT -1 65 "T i
s the maximum time to be computed; d is the diffusion constant" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "L:=1; k:=0.05; h:=0.10; T:=10; d:=0
.1;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 60 "r and c are the constants \+
appearing in the iteration formula" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
22 "r:= d*k/h^2; c:=1-2*r;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 10 "par
ameter " }{XPPEDIT 18 0 "lambda;" "6#%'lambdaG" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 10 "lambda:=5;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "No
nlinearity f(u)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "f:=y->eval(y*(1-
y)); plot(f(y), y=0..1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "Initi
al Condition" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "u0:=y->evalf(sin(Pi
*y)); plot(u0(y), y=0..L);" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 29 "Dirichlet boundary conditions" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 25 "left:=y->0; right:=y->0; " }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 26 "Generate the spatial grids" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 44 "n:=floor(L/h)+1; x:=[seq(0+(i-1)*h,i=1..n)];" }}}{EXCHG {PARA 
0 "" 0 "" {TEXT -1 24 "Generate the time grids " }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 44 "m:=floor(T/k)+1; t:=[seq(0+(j-1)*k,j=1..m)];" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 38 "Initialize the solution matrix by \+
zero" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "u:=matrix(n,m,[seq([seq(0.0
,j=1..m)],i=1..n)]):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 23 "Input the
 initial value" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "for i from 1 to n
 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "     u[i,1]:=u0(x[i]);" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 29 "Input the boundary conditions" }{MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for j from 1 to m-1 do" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "     u[1,j+1]:=left(t[j]); u[n,j+1]
:=right(t[j]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 43 "Compute the matrix by the iteration formu
la" }{MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "for j from
 1 to m-1 do  " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "   for i from 2 t
o n-1 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "      u[i, j+1]:=r*u[i-
1,j]+c*u[i,j]+r*u[i+1,j]+lambda*k*f(u[i,j]); " }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 11 "   end do; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "end \+
do:" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 9 "Animation" }}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 70 "Each of the following is a static plot of
 the solution at a given time" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 20 "for j from 1 to m do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "  U[j
]:=plot([seq([x[i],u[i,j]],i=1..n)]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 7 "end do:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 58 "display the seque
nce of the static plots into an animation" }{MPLTEXT 1 0 0 "" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 45 "display([seq(U[j],j=1..m)], insequence=tr
ue);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 
"" 0 "" {MPLTEXT 0 21 46 "**** THIS IS A PLOT THAT CAN BE ANIMATED ***
**" }}}}}{MARK "1 14 4 0" 3 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }
{PAGENUMBERS 0 1 2 33 1 1 }