#C20.txt, April 7, 2014
Help:=proc(): print(`Dirichlet(f,x,y,h,k,x1,y1);`): end:

#appx. value at (x1,y1) of
#the sol. of the elliptic bvp
#u_{xx}+u_{yy}=f(x,y) , in the unit square
#with Dirichlet boundary conditions (0 on the bondary)
#code written by Cole Franks
#(disclaimer: not repsonsible to any damage caused ...)
Dirichlet:=proc(f,x,y,h,k,x1,y1) local T,j,n,N,J,r,eq,var:

J:=1/h:
N:=1/k:

for j from 0 to J do
T[j,0]:=0: T[j,N]:=0:
od:

for n from 0 to N do
T[0,n]:=0: T[J,n]:=0:
od:

var:={seq(seq(T[j,n], j = 1..J-1), n = 1..N-1)}:

eq:={seq(seq(subs({x = h*j, y = k*n},f) =

(T[j+1,n]-2*T[j,n]+T[j-1,n])/h^2 + (T[j,n+1]-2*T[j,n]+T[j,n-1])/k^2,

j = 1..J-1), n = 1..N-1)}:

var:=solve(eq, var):


subs(var, T[x1/h, y1/k]):

#T[x1/h, y1/k]:



end: