function pde3b
% pde3b.m  3D non equal h boundary value PDE, iterative solution  
%  u(x,y,z) = x^3  + y^3  +z^3  + 2x^2y + 3xy^2 + 4z^2x + 5y + 6  
%     du/dx = 3x^2 + 0    + 0   + 4xy   + 3y^2  + 4z^2  + 0  + 0  
% d^2u/dx^2 = 6x   + 0    + 0   + 4y    + 0     + 0     + 0  + 0  
%     du/dy = 0    + 3y^2 + 0   + 2x^2  + 6xy   + 0     + 5  + 0  
% d^2u/dy^2 = 0    + 6y   + 0   + 0     + 6x    + 0     + 0  + 0  
%     du/dz = 0    + 0    +3z^2 + 0     + 0     + 8zx   + 0  + 0  
% d^2u/dz^2 = 0    + 0    +6z   + 0     + 0     + 8x    + 0  + 0  
%                                                        
% solve d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 = c(x,y,z)     
% by second order method on cube -1,-1,-1 to 1.2,1,1     

% second order numerical difference                      
% d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 =                    
% 1/hx^2(u(x-hx,y,z)+u(x+hx,y,z)-2u(x,y,z)) +            
% 1/hy^2(u(x,y-hy,z)+u(x,y+hy,z)-2u(x,y,z))              
% 1/hz^2(u(x,y,z-hz)+u(x,y,z+hz)-2u(x,y,z)) + O(h^2)     
%  = c(x,y,z) = 20x + 10y + 6z                           
% solve for new u(x,y,z) =                               
% ((u(x-hx,y,z)+u(x+hx,y,z))/hx^2 +                      
%  (u(x,y-hy,z)+u(x,y+hy,z))/hy^2 +                      
%  (u(x,y,z-hz)+u(x,y,z+hz))/hz^2 - c(x,y,z)) /          
%  (2/hx^2 + 2/hy^2 + 2/hz^2)                            

format compact; % prevent blank lines
fid = fopen('pde3b_m.out', 'w');
fprintf(fid, 'pde3b.m running \n');
sprintf('generating pde3b_m.out \n')
xmin = -1.0;
ymin = -1.0;
zmin = -1.0;
xmax =  1.2;
ymax =  1.0;
zmax =  1.0;
nx = 11;
ny = 11;
nz = 9;
hx   = (xmax-xmin)/(nx-1); 
hy   = (ymax-ymin)/(ny-1);
hz   = (zmax-zmin)/(nz-1);
us = zeros(nx, ny, nz); % initialize matrix
ut = zeros(nx, ny, nz);
diff = 1.0E100;     % declare variables then we
error = 0.0;        % do not need 'global'
% avgerror = 0.0;   % only used locally, lint complains
maxerror = 0.0;
iter = 0;

fprintf(fid, 'xmin=%g, xmax=%g, nx=%d, hx=%g \n', xmin, xmax, nx, hx);
fprintf(fid, 'ymin=%g, ymax=%g, ny=%d, hy=%g \n', ymin, ymax, ny, hy);
fprintf(fid, 'zmin=%g, zmax=%g, nz=%d, hz=%g \n', zmin, zmax, nz, hz);
init_boundary();
fprintf(fid, '\n initial conditions set \n');

while diff>0.0001 && iter<149
  iterate();
  check();
  avgerror = error/((nx-2).*(ny-2).*(nz-2));
  fprintf(fid, 'iter=%4d, diff=%8.4f, error=%9.4f, avg=%8.4f, max=%8.4f\n',...
          iter, diff, error, avgerror, maxerror);
  iter = iter+1;
  if mod(iter,50)==0
    printus();
    printdiff();
  end
end
printexact()
printus()
printdiff()
fclose(fid);
return;

function value = u(x,y,z)
  % solution for boundary and checking
  value = x.*x.*x + y.*y.*y + z.*z.*z + 2.0.*x.*x.*y + 3.0.*x.*y.*y + ...
          4.0.*z.*z.*x + 5.0.*y + 6.0;
end % z

function term = c(x, y, z)
  % the term on the right hand side
  term = 20.0.*x + 10.0.*y + 6.0.*z;
end % c

function value = u000(i, j, k)
  % compute new iterated value of us(i,j,k)
  % the equation is the solution of the difference solution
  value = ((us(i-1,j,k) + us(i+1,j,k))/(hx*hx) + ... 
           (us(i,j-1,k) + us(i,j+1,k))/(hy*hy) + ...
           (us(i,j,k-1) + us(i,j,k+1))/(hz*hz) - ...
           c(xmin+(i-1)*hx,ymin+(j-1)*hy,zmin+(k-1)*hz)) / ...
           (2.0/(hx*hx) + 2.0/(hy*hy) + 2.0/(hz*hz));
end % u000

function init_boundary()
  % set boundary on six sides 
  for i=1:nx
    for j=1:ny
      us(i,j,1)  = u(xmin+(i-1)*hx, ymin+(j-1)*hy, zmin);
      us(i,j,nz) = u(xmin+(i-1)*hx, ymin+(j-1)*hy, zmax);
    end
  end
  for j=1:ny
    for k=1:nz
      us(1,j,k)  = u(xmin, ymin+(j-1)*hy, zmin+(k-1)*hz);
      us(nx,j,k) = u(xmax, ymin+(j-1)*hy, zmin+(k-1)*hz);
    end
  end
  for i=1:nx
    for k=1:nz
      us(i,1,k)  = u(xmin+(i-1)*hx, ymin, zmin+(k-1)*hz);
      us(i,ny,k) = u(xmin+(i-1)*hx, ymax, zmin+(k-1)*hz);
    end
  end
end % init_boundary

function iterate()
  % perform one iteration step including copy
  diff = 0.0;
  for i=2:nx-1
    for j=2:ny-1
      for k=2:nz-1
        ut(i,j,k) = u000(i, j, k);
        diff = diff + abs(ut(i,j,k)-us(i,j,k));
      end
    end
  end
  % copy temp to solution
  for i=2:nx-1
    for j=2:ny-1
      for k=2:nz-1
        us(i,j,k) = ut(i,j,k);
      end
    end
  end
end % iterate

function check()
  % compute total error and maximum cell error
  error = 0.0;
  maxerror = 0.0;
  for i=2:nx-1
    for j=2:ny-1
      for k=2:nz-1
        zexact = u(xmin+(i-1)*hx, ymin+(j-1)*hy, zmin+(k-1)*hz);
        err = abs(us(i,j,k)-zexact);
        error = error + err;
        if err>maxerror
          maxerror = err;
        end
      end
    end
  end
end % check

function printus()
  % print the solution thus far
  fprintf(fid, '\n computed solution \n');
  for k = 2:nz-1
    fprintf(fid, 'u(x,y,%d) plane \n', k);
    for i=2:nx-1
      for j=2:ny-1
        fprintf(fid, '%7.3f', us(i,j,k));
      end
      fprintf(fid, '\n');
    end
    fprintf(fid, '\n');
  end
end % printus

function printexact()
  % print exact solution
  fprintf(fid, '\n exact solution \n');
  for k = 1:nz
    fprintf(fid, 'u(x,y,%d) plane \n', k);
    for i=1:nx
      for j=1:ny
        fprintf(fid, '%7.3f', u(xmin+(i-1)*hx, ymin+(j-1)*hy, zmin+(k-1)*hz));
      end
      fprintf(fid, '\n');
    end
    fprintf(fid, '\n');
  end
end % printexact

function printdiff()
  % print exact solution minus computed solution
  fprintf(fid, '\n difference exact-computed solution \n');
  for k=2:nz-1
    fprintf(fid, 'u(x,y,%d) plane \n', k);
    for i=2:nx-1
      fprintf(fid, '       ');
      for j=2:ny-1
        fprintf(fid, '%7.4f', u(xmin+(i-1)*hx, ymin+(j-1)*hy, ...
                     zmin+(k-1)*hz)-us(i,j,k));
      end
      fprintf(fid, '\n');
    end
    fprintf(fid, '\n');
  end
end % printdiff

end % pde3b.m