{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 56 "L:=1; k:=0.004; h:=0.10; T:=2; d1:=
0.00002; d2:=0.00001;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 108 "r and c
 are the constants appearing in the iteration formula   (make sure tha
t  r<0.5  for stable algorithm)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "
r1:= d1*k/h^2; r2:=d2*k/h^2; c1:=1-2*r1; c2:=1-2*r2;" }}}{EXCHG {PARA 
0 "" 0 "" {TEXT -1 30 "Nonlinearity f(u,v) and g(u,v)" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 76 "F:=0.03; K:=0.06; f:=(u,v)->0.03*(1-u)-u^2*v; g:
=(u,v)->u^2*v-(0.03+0.06)*v;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "I
nitial Condition" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 88 "u0:=y->evalf(0.
01*sin(Pi*y)); v0:=y->evalf(0.05*sin(Pi*y)); plot([u0(y),v0(y)], y=0..
L);" }{TEXT -1 1 " " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 26 "Generate t
he 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:=floo
r(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 95 "u:=matrix(n,m,[seq([seq(0.0,j=1..m)],i=1..n)]): v:
=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 40 "   \+
  u[i,1]:=u0(x[i]); v[i,1]:=v0(x[i]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 4 "end:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 43 "Compute the matrix b
y the iteration formula" }{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 109 "      u[1,j+1]:=2*r1*u[2,j]+c1*u[1,j]+k*f(u[1,j],v[1
,j]); v[1,j+1]:=2*r2*v[2,j]+c2*v[1,j]+k*g(u[1,j],v[1,j]);" }}{PARA 0 "
> " 0 "" {MPLTEXT 1 0 113 "      u[n,j+1]:=2*r1*u[n-1,j]+c1*u[n,j]+k*f
(u[n,j],v[n,j]); v[n,j+1]:=2*r2*v[n-1,j]+c2*v[n,j]+k*g(u[n,j],v[n,j]);
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "   for i from 2 to n-1 do" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 136 "      u[i, j+1]:=r1*u[i-1,j]+c1*u[
i,j]+r1*u[i+1,j]+k*f(u[i,j],v[i,j]); v[i, j+1]:=r2*v[i-1,j]+c2*v[i,j]+
r2*v[i+1,j]+k*g(u[i,j],v[i,j]); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 10 
"   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 giv
en time" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "for j from 1 to \+
m do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "  U[j]:=plot(\{[seq([x[i],u
[i,j]],i=1..n)], [seq([x[i],v[i,j]],i=1..n)]\});" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 7 "end do:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 58 "displa
y the sequence 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)], in
sequence=true);" }}}{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 "3" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 
}{PAGENUMBERS 0 1 2 33 1 1 }