function nma_P2DDIRJCB_S()
%script solves 2D Poission PDE on unit square using Jacobian iterative method

% file name: nma_P2DDIRJCB_S.m
%
% Small matlab script which solves the 2D Poission PDE
% on unit square using the Jacobian iterative method
% with stopping criteria that of the relative residual.
%
% This is a standalone script which require no additional
% toolboxes to run.
%
% There are 3 parameters to change only below. Edit these and rerun
% to see the solution
%
% n: THis is overall number of grid points on one edge.
% c: Tolerance constant. Use 1 or 0.1 or 0.01, used for checking for
%    convergence.
% force: this is the RHS of the PDE, edit to change the function
%
% by Nasser M. Abbasi  11/17/2010
%

%---------- CONFIGURATION PARAMETERS ---------------------------
n     = 9;   % all grid points, including internal points, on one edge
c     = 0.01; % constant of tolerance
myforce = @(X,Y) -exp( -(X-0.25).^2 - (Y-0.6).^2 ); % RHS of PDE
%---------- END CONFIGURATION PARAMETERS ------------------------

h     = 1/(n-1);        % grid spacing
tol   = c * h^2;        % tolerance
res   = zeros(n,n);     % storage for residual
u     = res;            % storage for solution
A     = [0 1 0;1  0 1;0 1 0];  % 5 point stencil for u calculation
B     = [0 1 0;1 -4 1;0 1 0];  % 5 point stencil for residue calculation
k     = 0;                     % loop counter
[X,Y] = meshgrid(0:h:1,0:h:1); % coordinates
f     = myforce(X,Y);
normf = norm( reshape(f(2:end-1,2:end-1),1,(n-2)^2),2); % find norm

% Now start the iteration solver, stop when relative residual < tolerance

figure;
done  = false;
i     = 2:n-1;
j     = 2:n-1;      %the iterations vectorized

while ~done
    k = k+1;
    
    u(i,j)   = (1/4)*( u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) - h^2 * f(i,j) );
    res(i,j) = f(i,j) - (1/h^2)*( u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) - 4*u(i,j) );
    
    mesh(X,Y,u);  hold on; title(sprintf('solution at iteration %d',k));
    %zlim([0,.07]);
    drawnow;
    
    if abs(normf)<eps
        rel = norm(reshape(res,1,n^2),2);
    else
        rel = norm(reshape(res,1,n^2),2)/normf;
    end
    
    if rel<tol
        done = true;
    end;
end
end