function pde3_eq
% pde3_eq.m  just checking second order approximation      
%      z(x,y) = x^3  + y^3  + 2x^2y + 3xy^2 + 4x + 5y + 6  
%       dz/dx = 3x^2 + 0    + 4xy   + 3y^2  + 4  + 0  + 0  
%   d^2z/dx^2 = 6x   + 0    + 4y    + 0     + 0  + 0  + 0  
%       dz/dy = 0    + 3y^2 + 2x^2  + 6xy   + 0  + 5  + 0  
%   d^2z/dy^2 = 0    + 6y   + 0     + 6x    + 0  + 0  + 0  
%                                                          
% solve d^2z/dx^2 + d^2z/dy^2 = 12x + 10y = c(x,y)         
% by second order method on square -1,-1 to 1,1            

% second order numerical difference                        
% d^2z/dx^2 + d^2z/dy^2 =                                  
% 1/h^2(z(-1h,0)+z(1h,0)+z(0,-1h)+z(0,1h)-4z(0,0)) +O(h^2) 
%  = c(x,y)                                                
% solve system of linear equations for                     
% z(i-1,j)+z(i+1,j)+z(i,j-1)+z(i,j+1)-4z(i,j) = c(x,y)*h^2 
% see pde3_eq.txt for development                          

format compact;
nx = 11; % conversion from "C" to Matlab 0..nx-1 to 1:nx
ny = 11; % the area inside the boundary is 2:nx-1
nn = (nx-2)*(ny-2);
zt = zeros(nn,nn); % simultaneous equations built here      
     % (nx-2)*(ny-2) by (nx-2)*(ny-2) internal, non boundary cells   
     % position is (i-1)+(nx-2)*(j-1) for internal solution z(i,j)   
     % (:,(nx-2)*(ny-2)) is constant,y  as in  A x = y, solve for x  
% zs = zeros(nn); % solution being computed
zc = zeros(nn); % constant terms in equation
xmin = -1.0;
ymin = -1.0;
xmax =  1.0;
ymax =  1.0;
h    =  (xmax-xmin)/(nx-1); % nx==ny only need one h 
error = 0.0;         % sum of absolute exact-computed         
%avgerror = 0.0;     % average error                          
maxerror = 0.0;      % maximum cell error                     

sprintf('output file pde3_eq_m.out being written \n')
fid = fopen('pde3_eq_m.out', 'w');
fprintf(fid, 'pde3_eq.m running\n'); 
fprintf(fid, 'xmin=%g, xmax=%g, nx=%d, h=%g \n', xmin, xmax, nx, h);
fprintf(fid, 'ymin=%g, ymax=%g, ny=%d, h=%g \n', ymin, ymax, ny, h);
  
init_matrix();
fprintf(fid, 'matrix initialized\n');
fprintf(fid, ' initial matrix, upper left \n');
for i=1:11
  for j=1:10;
    fprintf(fid, '%5.2f ',zt(i,j));
  end
  fprintf(fid, '=%5.2f \n', zc(i));
end 
fprintf(fid, '\n initial matrix, lower right \n');
for i=nn-11:nn
  for j=nn-10:nn
    fprintf(fid, '%5.2f ',zt(i,j));
  end
  fprintf(fid, '=%5.2f \n', zc(i));
end
fprintf(fid, '\n');
  
%  simeq((nx-2)*(ny-2), (nx-2)*(ny-2)+1, zt, zs);
zs = zt\zc; %  solves zt * zs = zc for zs
printexact();
printzs();
printdiff();
  
check();
avgerror = error/((nx-2)*(ny-2));
fprintf(fid, 'total error=%g, average error=%g, max error=%g\n', ...
             error, avgerror, maxerror);
return;

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

function value = c(x, y) % from solve 
  value = 12.0.*x + 10.0.*y;
end

function init_matrix() 
  fprintf(fid, 'internal cells zeroed, h=%g\n', h);          
  for i=2:nx-1 % inside boundary 
    for j=2:ny-1
      ii = (i-2)+(nx-2)*(j-2)+1; % row of matrix 
      fprintf(fid, 'i=%d, j=%d, ii=%d, ', i, j, ii);
      zt(ii,ii) = -4.0;
      zc(ii) = h*h*c(xmin+(i-1)*h, ymin+(j-1)*h);
      if (j+1)<ny
        jj = (i-2)+(nx-2)*((j-2)+1)+1;
        zt(ii,jj) = 1.0;
        fprintf(fid, 'jj+=%d, ', jj);
      else
         zc(ii) = zc(ii) - z(xmin+(i-1)*h, ymin+(j-1+1)*h);
      end
      if (j-1)>1
        jj = (i-2)+(nx-2)*((j-2)-1)+1;
        zt(ii,jj) = 1.0;
        fprintf(fid, 'jj-=%d, ', jj);
      else
         zc(ii) = zc(ii) - z(xmin+(i-1)*h, ymin+(j-1-1)*h);
      end
      if (i+1)<nx
         jj = ((i-2)+1)+(nx-2)*(j-2)+1;
         zt(ii,jj) = 1.0;
         fprintf(fid, 'jji+=%d, ', jj);
      else
         zc(ii) = zc(ii) - z(xmin+(i-1+1)*h, ymin+(j-1)*h);
      end
      if (i-1)>1
         jj = ((i-2)-1)+(nx-2)*(j-2)+1;
         zt(ii,jj) = 1.0;
         fprintf(fid, 'jji-=%d', jj);
      else
         zc(ii) = zc(ii) - z(xmin+(i-1-1)*h, ymin+(j-1)*h);
      end
      fprintf(fid, '\n');
    end
  end
end  % init_matrix

function check() % typically not known, instrumentation 
  error = 0.0;
  maxerror = 0.0;
  for i=2:nx-1
    for j=2:ny-1
      ii = (i-2)+(nx-2)*(j-2)+1;
      zexact = z(xmin+(i-1)*h,ymin+(j-1)*h);
      err = abs(zs(ii)-zexact);
      error = error + err;
      if err>maxerror
        maxerror=err;
      end
    end
  end
end  % check

function printzs()
  fprintf(fid, '\n computed solution\n');
  for i=2:nx-1
    for j=2:ny-1
      ii = (i-2)+(nx-2)*(j-2)+1;
      fprintf(fid, '%6.3f ', zs(ii));
    end
    fprintf(fid, '\n');
  end
  fprintf(fid, '\n');
end   % printzs

function printexact()
  fprintf(fid, '\n exact solution\n'); % inside boundary
  for i=2:nx-1
    for j=2:ny-1
      fprintf(fid, '%6.3f ', z(xmin+(i-1)*h,ymin+(j-1)*h));
    end
    fprintf(fid, '\n');
  end
  fprintf(fid, '\n');
end  % printexact

function printdiff()
  fprintf(fid, '\n difference exact-computed solution\n');
  fprintf(fid, '\n'); 
  for i=2:nx-1
    for j=2:ny-1
      ii = (i-2)+(nx-2)*(j-2)+1;
      fprintf(fid, '%7.4f ', z(xmin+(i-1)*h,ymin+(j-1)*h)-zs(ii));
    end
    fprintf(fid, '\n');
  end
  fprintf(fid, '\n');
end  % printdiff

end  % pde3_eq.m