% FIGURE 1.8.1
%                              JACOBI METHOD
      M = 10;
      H = 1.0/M;
%                              SET BOUNDARY KNOWNS TO ZERO PERMANENTLY
%                              AND INTERIOR UNKNOWNS TO ZERO TEMPORARILY
      for I=1:M+1
      for J=1:M+1
      for K=1:M+1
         UOLD(I,J,K) = 0.0;
      end
      end
      end
%                              BEGIN JACOBI ITERATION
      for ITER = 1:10000
%                              UPDATE UNKNOWNS ONLY
         for I=2:M
         for J=2:M
         for K=2:M
            UNEW(I,J,K)  = H^2/6.0 + ( UOLD(I+1,J,K) + UOLD(I-1,J,K)    ...
                                     + UOLD(I,J+1,K) + UOLD(I,J-1,K)    ...
                                     + UOLD(I,J,K+1) + UOLD(I,J,K-1))/6.0;
         end
         end
         end
%                              COPY UNEW ONTO UOLD
         for I=2:M
         for J=2:M
         for K=2:M
           UOLD(I,J,K) = UNEW(I,J,K);
         end
         end
         end
%                              EVERY 10 ITERATIONS CALCULATE MAXIMUM
%                              RESIDUAL AND CHECK FOR CONVERGENCE
         if (mod(ITER,10) ~= 0)
            continue
         end
         RMAX = 0.0;
         for I=2:M
         for J=2:M
         for K=2:M
            RESID = 6*UOLD(I,J,K) - UOLD(I+1,J,K) - UOLD(I-1,J,K)     ...
                                  - UOLD(I,J+1,K) - UOLD(I,J-1,K)     ...
                                  - UOLD(I,J,K+1) - UOLD(I,J,K-1) - H^2;
            RMAX = max(RMAX,abs(RESID));
         end
         end
         end
         ITER
         RMAX = RMAX/H^2
         if (RMAX <= 1.e-10)
            break
         end
      end