% BVP_2D_Poisson_FD.m
% This MATLAB routine solves a 2-D boundary value
% problem on a rectangular domain using finite differences.
% Zero Dirichlet boundary conditions are used.
% K. Beers. MIT ChE. 9/4/03
function iflag_main = ...
    BVP_2D_Poisson_FD(fun_name,L,H,num_pts);
iflag_main = 0;

% set number of points in grid
Nx = num_pts; Ny = num_pts;  Ntot = Nx*Ny;

% First, place a 2-D computational grid.
x = linspace(0,L,Nx); dx = x(2)-x(1);
y = linspace(0,H,Ny); dy = y(2)-y(1);

% Form a 2-D regular mesh such that
%     XG(i,j) = x(i), YG(i,j) = y(j)
[XG,YG] = form_2D_mesh(x,y);

% compute value of f(x,y) at each grid point
% and make a filled contour plot
FG = feval(fun_name,XG,YG,L,H);
figure; contourf(XG,YG,FG,min(Nx,Ny)); colorbar;
xlabel('x'); ylabel('y'); title('f(x,y)');

% allocate memory for the matrix and RHS vector
A = spalloc(Ntot,Ntot,5*Ntot); b = zeros(Ntot,1);

% We next specify equations for each boundary point.
% BC # 1
i = 1;
for j=1:Ny
    n = get_label(i,j,Nx,Ny);
    A(n,n) = 1;  b(n) = 0;
end
% BC # 2
i = Nx;
for j=1:Ny
    n = get_label(i,j,Nx,Ny);
    A(n,n) = 1;  b(n) = 0;
end
% BC # 3
j = 1;
for i=2:(Nx-1)
    n = get_label(i,j,Nx,Ny);
    A(n,n) = 1; b(n) = 0;
end
% BC # 4
j = Nx;
for i=2:(Nx-1)
    n = get_label(i,j,Nx,Ny);
    A(n,n) = 1; b(n) = 0;
end

% We now set the linear equations for the interior points.
factor_x = 1/(dx^2); factor_y = 1/(dy^2);
factor_cent = 2*(factor_x + factor_y);
for i=2:(Nx-1)
for j=2:(Ny-1)
    n = get_label(i,j,Nx,Ny);
    A(n,n-Ny) = -factor_x;  A(n,n+Ny) = -factor_x;
    A(n,n-1) = -factor_y;  A(n,n+1) = -factor_y;
    A(n,n) = factor_cent;
    b(n) = FG(i,j);
end
end

% We now solve using Gaussian elimination
phi = A\b;

% We now extract the results of the calculation
% into the 2-D grid format and make a filled
% contour plot.
PHIG = zeros(size(XG));
for i=1:Nx
for j=1:Ny
    n = get_label(i,j,Nx,Ny);
    PHIG(i,j) = phi(n);
end
end
figure; contourf(XG,YG,PHIG,min(Nx,Ny)); colorbar;
xlabel('x'); ylabel('y'); title('Poisson BVP (FD): \phi(x,y)');

% save the results to a .mat file
save BVP_2D_Poisson_FD.mat;

iflag_main = 1;
return;


%--------------------------------------
function [XG,YG] = form_2D_mesh(x,y);

Nx = length(x); Ny = length(y);
XG = zeros(Nx,Ny); YG = zeros(Nx,Ny);
for k=1:Nx
    XG(k,:) = x(k)*ones(1,Ny);
end
for k=1:Ny
    YG(:,k) = y(k)*ones(Nx,1);
end

return;


%--------------------------------------
function n = get_label(i,j,Nx,Ny);

n = (i-1).*Ny + j;

return;