%{
 Solving the Poisson equation fxx + fyy = F(x,y)
 Method:
 Finite-difference 2nd order
 AG: April 2022
To Do
# add derivative boundary conditions
%}

clear  % Clears command window
clc    % Clears command history
clf    % Removes anything in the figure window before simulation. 
tic    % set timer for the program

% INPUT
% Numerical control
iter = 10000;
tol = 1.0e-06;
omega = 1.2;
% Domain of the solutions
xmin = 0.0;
xmax = pi;
ymin = 0.0;
ymax = pi;
% the grid
Nx = 51;
Ny = 51;
dx = (xmax-xmin)/(Nx-1);
dy = (ymax-ymin)/(Ny-1);
x = zeros(1,Nx);
y = zeros(1,Ny);
f = zeros(Ny,Nx);
F = zeros(Ny,Nx);

% BOUNDARY conditions
% along x
for i = 1:Nx
  x(i) = xmin + dx*(i-1);
  f(i,1) = 1.0*sin(x(i));
  f(i,Ny)= abs(sin(2.0*x(i)));
end
% along y
for j = 1:Ny
  y(j) = ymin + dy*(j-1);
  f(1,j) = 0.0;
  f(Nx,j) = 0.0;
end
% Non-homogeneous term (Poisson equation)
F(ceil(Nx/2),ceil(Ny/2)) = 350.0;

% initial guess
% for i = 2:Nx-1
%   for j = 2:Ny-1
%     f(i,j) = 0.3*rand;
%   end
% end

figure (1)
surf(y,x,f)

% call the solver
[f2] = pde01(f,F,dx,dy,Nx,Ny,iter,tol,omega);
%===

figure (2)
surf(y,x,f2)
xlabel('y'), ylabel('x'), zlabel('Solution')

%print time of calculations
timeE = floor(toc);
fprintf('\n timeElapsed = %4i sec \n',timeE)

%{
  Solving the Poisson equation with Dirichlet BCs
  INPUT:
  f(i,j)    boundary conditions + initial guess for internal points
  F(i,j)    nonhomogeneous term
  dx, dy    grid increments
  Nx        number of grid points in x direction 
  Ny        number of grid points in y direction
  iter      maximum number of iterations
  tol       convergence tolerance
  omega     over-relaxation factor
  OUTPUT
  f2(x,y)   the solution
AG: April 2022
%}
function[f2] = pde01(f,F,dx,dy,Nx,Ny,iter,tol,omega)
% preparation
f2 = zeros(Nx,Ny);
beta2 = (dx/dy)^2;
d = 2.0*(1.0+beta2);
% iterations
for it=1:iter
  dfmax= 0.0;
  for j=2:Ny-1
    for i=2:Nx-1
      df=(f(i+1,j)+beta2*f(i,j+1)+f(i-1,j)+beta2*f(i,j-1)+dx^2*F(i,j)-d*f(i,j))/d; 
        if (abs(df) >= dfmax) 
          dfmax = df;
        end
      f(i,j)=f(i,j)+omega*df;
    end
  end
  if(dfmax <= tol) 
    fprintf(' \n The solution has converged, it = %3i',it)
    break
  end
end
% Prepare OUTPUT 
for j=1:Ny
  for i=1:Nx
    f2(i,j)= f(i,j);
  end
end

if it==iter
  fprintf('  \n The solution failed to converge, iter = %3i',iter)
end

end