% FIGURE 1.8.2
%                              GAUSS-SEIDEL 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
         U(I,J,K) = 0.0;
      end
      end
      end
%                              BEGIN GAUSS-SEIDEL ITERATION
      for ITER = 1:10000
%                              UPDATE UNKNOWNS ONLY
         for I=2:M
         for J=2:M
         for K=2:M
            GAUSS = H^2/6.0 + ( U(I+1,J,K) + U(I-1,J,K)     ...
                              + U(I,J+1,K) + U(I,J-1,K)     ...
                              + U(I,J,K+1) + U(I,J,K-1))/6.0;
            U(I,J,K) = GAUSS;
         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*U(I,J,K) - U(I+1,J,K) - U(I-1,J,K)     ...
                               - U(I,J+1,K) - U(I,J-1,K)     ...
                               - U(I,J,K+1) - U(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