%Who  : Abdallah Sayyed Ahmad.
%What : 2D Poission equation solver;DelV=f and bc1 at the bottom, bc2 at the right, bc3 at top and bc4 at left of a rectangular box
%Where: Computational Physics 338/Department of Physics/Birzeit University
%When : December 2020

clear;
clc;
f = inline('-2*((1-6*x^2)*y^2*(1-y^2)+(1-6*y^2)*x^2*(1-x^2))', 'x', 'y');
bc1 = inline('0', 'x');
bc2 = inline('0', 'y');
bc3 = inline('0', 'x');
bc4 = inline('0', 'y');

N=100;
x0=0;
x1=10;
y0=x0;
y1=x1;
dx=(x1-x0)/(N-1);
dy=dx;
b=zeros(N*N,1);
A=-1*gallery("poisson",N);


for j=2:N-1
      yy=y0+(j-1)*dy;      
  for i=2:N-1
    xx=x0+(i-1)*dx;
    indx=i+(j-1)*N;
    b(indx)=f(xx,yy)*dx*dx;
  end
end

% boundary conditions
for i=1:N
  %bottom
  indx=i;
  xx=x0+(i-1)*dx;
  b(indx)=bc1(xx);
  A(indx,:)=0;
  A(indx,indx)=1;

  %top
  indx=i+N*(N-1);
  xx=x0+(i-1)*dx;
  b(indx)=bc3(xx);
  A(indx,:)=0;
  A(indx,indx)=1;
end

for j=1:N
 %left
 indx=1+(j-1)*N;
  yy=y0+(j-1)*dy;
  b(indx)=bc4(yy);
  A(indx,:)=0;
   A(indx,indx)=1;

 %right
 indx=N+(j-1)*N;
  yy=y0+(j-1)*dy;
  b(indx)=bc2(yy);
  A(indx,:)=0;
   A(indx,indx)=1;

end
tic
sol=A\b;
%sol= bicgstab(A,b,1e-6,1000); 
toc
x = linspace(x0,x1,N);
y = linspace(y0,y1,N);
[xx,yy] = meshgrid(x,y);
Z = reshape(sol(:), N, []).';
%draw solution
contourf(xx,yy,Z)
% compute exact solution; this need to be changed when the source term or the boundary conditions changed
%exact=yy.*(1-yy).*xx.^3;
% compute error in solution per grid point
%error=norm(Z-exact)/N^2