function outvar = Poisson_2D(D,S,BC,a,b,c,d,m,n)
% Poisson_2D: 2d Poisson solver from Frederic Gibou

tic; % To compute the execution time.

% Discretize the domain:
x = linspace(a,b,m);
y = linspace(c,d,n);
dx = (b-a)/(m-1);
dy = (d-c)/(n-1);
mn = m*n;

% Initialize the variables:

A = sparse(mn,mn); % This is more computationally efficient than A = zeros(mn,mn).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assemble the linear system: %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Treat the interior points:

for i = 2:m-1
    for j = 2:n-1
        r = (j-1)*m+i;
        A(r,r  ) =  2*D/dx/dx + 2*D/dy/dy;
        A(r,r-1) =  -D/dx/dx;
        A(r,r+1) =  -D/dx/dx;
        A(r,r-m) =  -D/dy/dy;
        A(r,r+m) =  -D/dy/dy;
        rhs(r) = S(x(i),y(j));
    end
end

% Treat the boundary conditions:

for j = 1:n
    i = 1; % Left boundary condition
    r = (j-1)*m+i;
    A(r,r) = 1; rhs(r) = BC(a,y(j));
    i = m; % Right boundary condition
    r = (j-1)*m+i;
    A(r,r) = 1; rhs(r) = BC(b,y(j));
end
for i = 1:m
    j = 1; % Bottom boundary condition
    r = (j-1)*m+i;
    A(r,r) = 1; rhs(r) = BC(x(i),c);
    j = n; % Top boundary condition
    r = (j-1)*m+i;
    A(r,r) = 1; rhs(r) = BC(x(i),d);
end

% Solve the linear system:
U = A\rhs';

toc; % To compute the execution time.

% Plot the solution:
for i = 1:m
    for j = 1:n
        exact(i,j) = BC(x(i),y(j));
    end
end
subplot(2,1,1);
surfl(x,y,exact); shading interp; colormap(summer);
axis([a b, c d, min(min(exact)) max(max(exact))]);
xlabel('x');ylabel('y'); zlabel('Temperature');
title('Exact Solution of the Temperature distribution');
subplot(2,1,2);

% Need to put Un, which is now a vector, into a grid form.
for i = 1:m
    for j = 1:n
        r = (j-1)*m+i;
        U_grid(i,j) = U(r);
    end
end
surfl(x,y,U_grid); shading interp; colormap(summer);
axis([a b, c d, min(min(exact)) max(max(exact))]);
xlabel('x');ylabel('y'); zlabel('Temperature');
title('Numerical Solution of the Temperature distribution');
outvar = U;
end
