% FIGURE 4.8.1
%                              N = NUMBER OF POINTS IN X,Y DIRECTIONS
%                              A = PDE PARAMETER
      N = 40;
      A = 40.0;
      DX = 1.0/N;
      DY = DX;
%                              U = 0
      for I=2:N
      for J=2:N
         U(I,J) = 0.0;
      end
      end
%                              SET VALUES OF U,P ON BOUNDARY (THESE WILL
%                              NEVER CHANGE)
      for L=1:N+1
         X(L) = (L-1)*DX;
         Y(L) = (L-1)*DY;
         U(L,1) = 0.0;
         P(L,1) = 0.0;
         U(L,N+1) = X(L);
         P(L,N+1) = 0.0;
         U(1,L) = 0.0;
         P(1,L) = 0.0;
         U(N+1,L) = Y(L)^3;
         P(N+1,L) = 0.0;
      end
%                              R = B - A*U
%                              P = R
      for I=2:N
      for J=2:N
         R(I,J) = (U(I+1,J)-2*U(I,J)+U(I-1,J))/DX^2 ...
                 +(U(I,J+1)-2*U(I,J)+U(I,J-1))/DY^2 ...
                 - A*U(I,J) + X(I)*Y(J)*(A*Y(J)^2-6);
         P(I,J) = R(I,J);
      end
      end
%                              BEGIN CONJUGATE GRADIENT ITERATION
      for ITER=1:1000
%                              AP = A*P
         for I=2:N
         for J=2:N
            AP(I,J) = (-P(I+1,J)+2*P(I,J)-P(I-1,J))/DX^2  ...
                     +(-P(I,J+1)+2*P(I,J)-P(I,J-1))/DY^2  ...
                     + A*P(I,J);
         end
         end
%                              LAMBDA = R*P / P*AP
         LAMBDA = DOT(R,P,N)/DOT(P,AP,N);
%                              U = U + LAMBDA*P
         U = UPDATE(U,U,LAMBDA,P,N);
%                              R = R - LAMBDA*AP
         R = UPDATE(R,R,-LAMBDA,AP,N);
%                              ALPHA = -R*AP / P*AP
         ALPHA = -DOT(R,AP,N)/DOT(P,AP,N);
%                              P = R + ALPHA*P
         P = UPDATE(P,R,ALPHA,P,N);
%                              CALCULATE MAXIMUM ERROR
         ERMAX = 0.0;
         for I=2:N
         for J=2:N
            ERR = abs(U(I,J) - X(I)*Y(J)^3);
            ERMAX = max(ERMAX,ERR);
         end
         end
         MAXITR = ITER;
%                              STOP IF ERROR < 10^(-6)
         if (ERMAX < 1.E-6)
            break
         end
      end
      MAXITR