function pde3 %  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(x-h,y)-2z(x,y)+z(x+h,y))+
% 1/h^2(z(x,y-h)-2z(x,y)+z(x,y+h)) + O(h^2)
%  = c(x,y)
% solve for new z(x,y) =
% (z(x-h,y)+z(x+h,y)+z(x,y-h)+z(x,y+h)-c(x,y)*h^2)/4
format compact; % prevent blank lines
fid = fopen('pde3_m.out', 'w');
fprintf(fid, 'pde3.m running \n');
sprintf('generating pde3_m.out \n')
xmin = -1.0;
xmax = 1.0;
ymin = -1.0;
ymax = 1.0;
nx = 11;
ny = 11;
h = (xmax-xmin)/(nx-1); % also = (ymax-ymin)/(ny-1)
zs = zeros(nx, ny); % initialize matrix
zt = zeros(nx, ny);
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, h=%g \n', xmin, xmax, nx, h);
fprintf(fid, 'ymin=%g, ymax=%g, ny=%d, h=%g \n', ymin, ymax, ny, h);
printexact();
init_boundary();
fprintf(fid, '\n initial condition \n');
printzs();

while diff>0.0001 && iter<200
  iterate();
  check();
  avgerror = error/((nx-2).*(ny-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
    printzs();
    printdiff();
  end
end
fclose(fid);
return;

function val = z(x,y)
  % solution for boundary and checking
  val = 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 % z

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

function value = z00(i, j)
  % compute new iterated value of zs(i,j)
  % the equation is the solution of the difference solution
  value = zs(i-1,j) + zs(i+1,j) + zs(i,j-1) + zs(i,j+1);
  value = (value - c(xmin+(i-1)*h, ymin+(j-1)*h).*h.*h)/4.0;
end % z00

function init_boundary()
  % initialize the boundary cells in zs(i,j)
  for i=1:nx
    zs(i,1)  = z(xmin+(i-1)*h, ymin);
    zs(i,ny) = z(xmin+(i-1)*h, ymax);
  end
  for j=1:ny
    zs(1,j)  = z(xmin, ymin+(j-1)*h);
    zs(nx,j) = z(xmax, ymin+(j-1)*h);
  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
      zt(i,j) = z00(i, j);
      diff = diff + abs(zt(i,j)-zs(i,j));
    end
  end
  % copy temp to solution
  for i=2:nx-1
    for j=2:ny-1
      zs(i,j) = zt(i,j);
    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
      zexact = z(xmin+(i-1)*h, ymin+(j-1)*h);
      err = abs(zs(i,j)-zexact);
      error = error + err;
      if err>maxerror
        maxerror = err;
      end
    end
  end
end % check

function printzs()
  % print the solution thus far
  fprintf(fid, '\n computed solution \n');
  for i=1:nx
    for j=1:ny
      fprintf(fid, '%7.3f', zs(i,j));
    end
    fprintf(fid, '\n');
  end
  fprintf(fid, '\n');
end % printzs

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

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

end % pde3.m