%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('0', 'x', 'y');
bc1 = inline('1', 'x');
bc2 = inline('0', 'y');
bc3 = inline('1', 'x');
bc4 = inline('1', 'y');

%k=1;

ctr=1;
for i=0:0.01:10
k_array(ctr)=i;
ctr=ctr+1;
end

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 ii=1:ctr-1
k=k_array(ii);
A(A==-4)=-4-((k^2)*(dx^2));


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, []).';


resV(ii)=Z(5,5);
end
%draw solution

%contourf(xx,yy,Z)
%surf(xx,yy,Z);

plot(k_array,resV), xlabel('k'), ylabel('V(5,5)');

% 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








