(* POISSON EQUATION DIFFERENCE ALGORITHM 12.1
*
*  To approximate the solution to the poisson equation
*    DEL(u) = F(x,y), a <= x <= b, c <= y <= d,
*    Subject to boundary conditions:
*           u(x,y) = G(x,y),
*	if x = a or x = b for c <= y <= d,
*	if y = c or y = d for a <= x <= b
*
*   INPUT: endpoints a, b, c, d; integers m, n; tolerance TOL;
*     	    maximum number of iterations NN
*
*  OUTPUT: Approximations W(i,J) to u(X(i),Y(J)) for each
*	    i = 1, ..., n-1 and J = 1, ..., M-1 or a message that 
*	    the maximum number of iterations was exceeded.
*)
TEMP = Input["This is the Finite-Difference method for Elliptic\n
  equations.\n
  Input the function F(x,y) in terms of x and y\n
    For example: x*Exp[y]\n"];
F[x_,y_] := Evaluate[TEMP];
TEMP=Input["Input the function G(x,y) in terms of x and y\n
    For example: x*Exp[y]\n"];
G[x_,y_] := Evaluate[TEMP];
OK = 0;
While[OK == 0,
   A = Input["Input the endpoint A of interval [A,B] on X-axis\n"];
   B = Input["Input the endpoint B of interval [A,B] on X-axis\n"];
   c = Input["Input the endpoint C of interval [C,D] on Y-axis\n"];
   d = Input["Input the endpoint D of interval [C,D] on Y-axis\n"];
   If[A >= B || c >= d,
      Input["Left endpoint must be less than right endpoint\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
OK = 0;
While[OK == 0,
   n = Input["Input the number or intervals N on the X-axis\n
   Note that N should be larger than 2\n"];  
   M = Input["Input the number or intervals M on the Y-axis\n
   Note that M should be larger than 2\n"];
   If[M <= 2 || n <= 2,
      Input["Numbers must exceed two\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
OK = 0;
While[OK == 0,
   TOL = Input["Input the tolerance.\n"];
   If[TOL > 0,
      OK = 1,
      Input["Tolerance must be positive\n
      \n
      Press 1 [enter] to continue\n"];
   ];
];
OK = 0;
While[OK == 0,
   NN = Input["Input maximum number of iterations\n
     no decimal points\n"];
   If[NN <= 0,
      Input["Must be a positive integer.\n
      \n
      Enter 1 [enter] to continue\n"],
      OK = 1;
   ];
];
If[OK == 1,
   M1 = M-1;
   M2 = M-2;
   N1 = n-1;
   N2 = n-2;
   (* Step 1 *)
   H = N[(B-A)/n];
   K = N[(d-c)/M];
   (* Steps 2 & 3 construct mesh points *)
   (* Step 2 *)
   For[i = 0, i <= n, i++,
      X[i] = N[A+i*H];
   ];
   (* Step 3 *)
   For[J = 0, J <= M, J++,
      Y[J] = N[c+J*K];
   ];
   (* Step 4 *)
   For[i = 1, i <= N1, i++,
      W[i,0] = N[G[X[i],c]];
      W[i,M] = N[G[X[i],d]];
   ];
   For[J = 0, J <= M, J++,
      W[0,J] = N[G[A,Y[J]]];
      W[n,J] = N[G[B,Y[J]]];
   ];
   For[i = 1, i <= N1, i++,
      For[J = 1, J <= M1, J++,
	 W[i,J] = 0;
      ];
   ];
   (* Step 5 - Use V for lambda, VV for mu *)
   V = N[H^2/K^2];
   VV = N[2*(1+V)];
   L = 1;
   OK = 0;
   (* Z is a new value of W(I,J) to be used in computing 
      the norm of the error E used in place of NORM *)
   (* Step 6 *)
   While[L <= NN && OK == 0,
      (* Steps 7-20 perform Gauss-Seidel iterations *)
      (* Step 7 *)
      Z = (-H^2*N[F[X[1],Y[M1]]]+N[G[A,Y[M1]]]+
           V*N[G[X[1],d]]+V*W[1,M2]+W[2,M1])/VV;
      e = Abs[W[1,M1]-Z];
      W[1,M1] = Z;
      (* Step 8 *)
      For[i = 2, i <= N2, i++,
         Z = N[(-H^2*F[X[i],Y[M1]]+V*G[X[i],d]+
             W[i-1,M1]+W[i+1,M1]+V*W[i,M2])/VV];
	 If[Abs[W[i,M1]-Z] > e,
	    e = Abs[W[i,M1]-Z];
         ];
	 W[i,M1] = Z;
      ];
      (* Step 9 *)
      Z = N[(-H^2*F[X[N1],Y[M1]]+G[B,Y[M1]]+
           V*G[X[N1],d]+W[N2,M1]+V*W[N1,M2])/VV];
      If[Abs[W[N1,M1]-Z] > e,
	 e = Abs[W[N1,M1]-Z];
      ];
      W[N1,M1] = Z;
      (* Step 10 *)
      For[LL = 2, LL <= M2, LL++,
	 J = M2-LL+2;
         (* Step 11 *)
	 Z = N[(-H^2*F[X[1],Y[J]]+G[A,Y[J]]+
               V*W[1,J+1]+V*W[1,J-1]+W[2,J])/VV];
         If[Abs[W[1,J]-Z] > e,
 	    e = Abs[W[1,J]-Z];
	 ];
	 W[1,J] = Z;
         (* Step 12 *)
         For[i = 2, i <= N2, i++,
   	    Z = N[(-H^2*F[X[i],Y[J]]+W[i-1,J]+
               V*W[i,J+1]+V*W[i,J-1]+W[i+1,J])/VV];
            If[Abs[W[i,J]-Z] > e,
	       e = Abs[W[i,J]-Z];
            ];
            W[i,J] = Z;
         ];
         (* Step 13 *)
         Z = N[(-H^2*F[X[N1],Y[J]]+G[B,Y[J]]+
              W[N2,J]+V*W[N1,J+1]+V*W[N1,J-1])/VV];
         If[Abs[W[N1,J]-Z] > e,
	    e = Abs[W[N1,J]-Z];
         ];
         W[N1,J] = Z;
      ];
      (* Step 14 *)
      Z = N[(-H^2*F[X[1],Y[1]]+V*G[X[1],c]+
            G[A,Y[1]]+V*W[1,2]+W[2,1])/VV];
      If[Abs[W[1,1]-Z] > e,
         e = Abs[W[1,1]-Z];
      ];
      W[1,1] = Z;
      (* Step 15 *)
      For[i = 2, i <= N2, i++,
         Z = N[(-H^2*F[X[i],Y[1]]+V*G[X[i],c]
             +W[i+1,1]+W[i-1,1]+V*W[i,2])/VV];
         If[Abs[W[i,1]-Z] > e,
	    e = Abs[W[i,1]-Z];
         ];
         W[i,1] = Z;
      ];
      (* Step 16 *)
      Z = N[(-H^2*F[X[N1],Y[1]]+V*G[X[N1],c]+
            G[B,Y[1]]+W[N2,1]+V*W[N1,2])/VV];
      If[Abs[W[N1,1]-Z] > e,
         e = Abs[W[N1,1]-Z];
      ];
      W[N1,1] = Z;
      (* Step 17 *)
      If[e <= TOL,
         (* Step 18 *)
         FLAG = Input["Select output destination\n
         1. Screen\n
         2. Text file\n
         Enter 1 or 2\n"];
         If[FLAG == 2,
            NAME = InputString["Input the file name\n
            For example:   output.dta\n"];
            OUP =
             OpenWrite[NAME,FormatType->OutputForm],
            OUP = "stdout";
         ];
         Write[OUP,"POISSON EQUATION FINITE-DIFFERENCE METHOD\n"];
         Write[OUP,"\n"];
         Write[OUP,"i   J   X(i)    Y(J)      W(i,J)\n"];
         Write[OUP,"\n"];
         For[i = 1, i <= N1, i++,
            For[J = 1, J <= M1, J++,
   	       Write[OUP,i,"  ",J,"    ",SetPrecision[X[i],10],
                      "    ",SetPrecision[Y[J],10],"        ",
                           SetPrecision[W[i,J],10]];
            ];
         ];
         Write[OUP,"Convergence occured on iteration number : \n",L];
         (* Step 19 *)
         OK = 1,
         (* Step 20 *)
         L = L+1;
      ];
   ];
   (* Step 21 *)
   If[L > NN,
      Print["Method fails after iteration number \n",NN];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP];
   ];
];
Quit[];