function pde3b_eq
% pde3b_eq.c  3D non equal h boundary value PDE  non iterative    
%  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) = 12x + 10y - 8z                                   

% solve system of linear equations for                           
% u(i-1,j,k)/hx^2 + u(i+1,j,k)/hx^2 +                            
% u(i,j-1,k)/hy^2 + u(i,j+1,k)/hy^2 +                            
% u(i,j,k-1)/hz^2 + u(i,j,k+1)/hz^2 -                            
% (2/hx^2 + 2/hy^2 + 2/hz^2)*u(i,j,k) = c(x,y,z)                 
% for i=1,nx-2 j=1,ny-2 k=1,nz-2 while substituting              
% u(i,j,k) boundary values i=0 or nx-1, j=0 or ny-1, k=0 or nz-1 
% x=xmin+i*hx, y=ymin+j*hz, z=zmin+k*hz                          

% create two dimensional array with                              
% (nx-2)*(ny-2)*(nz-2) columns, one for each equation            
% and (nx-2)*(ny-2)*(nz-2)+1 rows for variables plus constant    
% Initialize the matrix, see init_matrix in the source code      
% Then solve the system of linear equations to get the solution. 

format compact;
nx = 11;
ny = 11;
nz = 9;
nn = (nx-2)*(ny-2)*(nz-2); % number of equations for non boundary
us = zeros(nn); % solution being computed 9*9*7      
ut = zeros(nn,nn); % temporary equations           
                   % nx-2 by ny-2 by nz-2 internal, non boundary cells
uc = zeros(nn);    % = constant in equation  us = ut\uc  or  ut*us=uc
xmin = -1.0;
ymin = -1.0;
zmin = -1.0;
xmax =  1.2;
ymax =  1.0;
zmax =  1.0;
hx   = (xmax-xmin)/(nx-1); 
hy   = (ymax-ymin)/(ny-1); 
hz   = (zmax-zmin)/(nz-1); 
error = 0.0;         % sum of absolute exact-computed         
%avgerror = 0.0;      % average error                          
maxerror = 0.0;      % maximum cell error                     

  sprintf('results in file  pde3b_eq_m.out')
  fid = fopen('pde3b_eq_m.out', 'w');  
  fprintf(fid, 'pde3b_eq.c running\n');
  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_matrix();
  fprintf(fid, 'matrix initialized\n');
  fprintf(fid, ' initial matrix, left side, upper \n');
  for i=1:11
    for j=1:11
      fprintf(fid, '%5.2f ',ut(i,j));
    end
    fprintf(fid, '\n');
  end
  fprintf(fid, '\n initial matrix, right side, upper \n');
  for i=1:11
    for j=nn-9:nn
      fprintf(fid, '%5.2f ',ut(i,j));
    end
    fprintf(fid, '=%5.2f\n', uc(i));
  end
  fprintf(fid, '\n');
  fprintf(fid, '\n initial matrix, right side, lower \n');
  for i=nn-10:nn
    for j=nn-9:nn
      fprintf(fid, '%5.2f ',ut(i,j));
    end
    fprintf(fid, '=%5.2f\n', uc(i));
  end
  fprintf(fid, '\n');
  
  %simeq((nx-2)*(ny-2)*(nz-2), (nx-2)*(ny-2)*(nz-2)+1, ut, us);
  us = ut\uc;
  printexact();
  printus();
  printdiff();

  check();
  avgerror = error/((nx-2)*(ny-2)*(nz-2));
  fprintf(fid, 'total error=%g, average error=%g, max error=%g\n', ...
         error, avgerror, maxerror);
  return
  
function value = u(x, y, z) % for boundary and check 
  % solution is not typically known
  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

function value = c(x, y, z) % RHS from solve 
  value = 20.0.*x + 10.0.*y + 6.0.*z;
end

function check() % typically not known, instrumentation 
  error = 0.0;
  maxerror = 0.0;
  for i=2:nx-1
    for j=2:ny-1
      for k=2:nz-1
        uexact = u(xmin+(i-1)*hx, ymin+(j-1)*hy, zmin+(k-1)*hz);
        ii = (i-2)+(nx-2)*(j-2)+(nx-2)*(ny-2)*(k-2)+1;
        err = abs(us(ii)-uexact);
        error = error + err;
        if err>maxerror
          maxerror=err;
        end
      end
    end
  end
end  % check 

function printus() 
  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
        ii = (i-2)+(nx-2)*(j-2)+(nx-2)*(ny-2)*(k-2)+1;
        fprintf(fid, '%6.3f ', us(ii));
      end
      fprintf(fid, '\n');
    end
  end
  fprintf(fid, '\n');
end  % printus 

function printexact() % typically not known 
  fprintf(fid, '\n exact 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, '%6.3f ', u(xmin+(i-1)*hx, ymin+(j-1)*hy, zmin+(k-1)*hz));
      end
      fprintf(fid, '\n');
    end
  end
  fprintf(fid, '\n');
end  % printexact 

function printdiff()
  fprintf(fid, '\n difference exact-computed solution\n');
  fprintf(fid, '\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
        ii = (i-2)+(nx-2)*(j-2)+(nx-2)*(ny-2)*(k-2)+1;
        fprintf(fid, '%7.4f', u(xmin+(i-1)*hx, ymin+(j-1)*hy, zmin+(k-1)*hz)-us(ii));
      end
      fprintf(fid, '\n');
    end
    fprintf(fid, '\n');
  end
end  % printdiff 

function init_matrix()
  hx2 = hx*hx;
  hy2 = hy*hy;
  hz2 = hz*hz; 
  
  fprintf(fid, 'internal cells initialized to zero \n');          
  for i=2:nx-1 % inside boundary 
    for j=2:ny-1
      for k=2:nz-1
        ii = (i-2)+(nx-2)*(j-2)+(nx-2)*(ny-2)*(k-2)+1; % column of matrix 
        ut(ii,ii) = -(2.0./hx2 + 2.0./hy2 + 2.0./hz2);
        uc(ii) = c(xmin+(i-1)*hx, ymin+(j-1)*hy, zmin+(k-1)*hz);
        if (i+1)<(nx-1)
          jj = ((i-2)+1)+(nx-2)*(j-2)+(nx-2)*(ny-2)*(k-2)+1;
          ut(ii,jj) = 1.0/hx2;
        else
          uc(ii) = uc(ii) - u(xmin+(i-1+1)*hx, ymin+(j-1)*hy, zmin+(k-1)*hz)/hx2;
        end
        if (i-1)>1
          jj = ((i-2)-1)+(nx-2)*(j-2)+(nx-2)*(ny-2)*(k-2)+1;
          ut(ii,jj) = 1.0/hx2;
        else
          uc(ii) = uc(ii) - u(xmin+(i-1-1)*hx, ymin+(j-1)*hy, zmin+(k-1)*hz)/hx2;
        end
        if (j+1)<(ny-1)
          jj = (i-2)+(nx-2)*((j-2)+1)+(nx-2)*(ny-2)*(k-2)+1;
          ut(ii,jj) = 1.0/hy2;
        else
          uc(ii) = uc(ii) - u(xmin+(i-1)*hx, ymin+((j-1)+1)*hy, zmin+(k-1)*hz)/hy2;
        end
        if (j-1)>1
          jj = (i-2)+(nx-2)*((j-2)-1)+(nx-2)*(ny-2)*(k-2)+1;
          ut(ii,jj) = 1.0/hy2;
        else
          uc(ii) = uc(ii) - u(xmin+(i-1)*hx, ymin+((j-1)-1)*hy, zmin+(k-1)*hz)/hy2;
        end
        if (k+1)<(nz-1)
          jj = (i-2)+(nx-2)*(j-2)+(nx-2)*(ny-2)*((k-2)+1)+1;
          ut(ii,jj) = 1.0/hz2;
        else
          uc(ii) = uc(ii) - u(xmin+(i-1)*hx, ymin+(j-1)*hy, zmin+((k-1)+1)*hz)/hz2;
        end
        if (k-1)>1
          jj = (i-2)+(nx-2)*(j-2)+(nx-2)*(ny-2)*((k-2)-1)+1;
          ut(ii,jj) = 1.0/hz2;
        else
          uc(ii) = uc(ii) - u(xmin+(i-1)*hx, ymin+(j-1)*hy, zmin+((k-1)-1)*hz)/hz2;
        end
      end
    end
  end
end   % init_matrix 
         
end  % pde3b_eq.m