%%
%Laplace Equation BVP example solution by iteration for 
%Cheng Example 4-7
%EP440 Example #2 of 2 for Numerical Project #1

%%
%Set domain
%Clear variables first
clear;
%Numerical Domain Size
Nx=100;
Ny=100;
%Physical Domain Size
a=1;
b=1;
%Resulting Domain Resolution
dx=a/(Nx-1);
dy=b/(Ny-1);
%Scale References
x=0:dx:a;
y=0:dy:b;

%%
%Construct numerical solution with Applied Boundary Conditions
%Number of Iterations
Niterations=10000;
%i and j domains (x and y assignment ranges)
i=2:Nx-1;
j=2:Ny-1;
%Initialize solution array
V=zeros(Nx,Ny);
%Specify boundary conditions – These will not be changed by the solver.
    V(:,1)=0;
    V(1,:)=1;
    V(:,Ny)=0;
    V(Nx,:)=0;
%Iterate Solution
for n=1:Niterations
    Vp=V;
    %Effectively force solution at any point to be a spatial average of
    %vertically and horizontally adjacent points, weighted by mesh scale.
    %By enforcing BCs at each step, we guarantee (eventual) convergence
    V(i,j)=((Vp(i+1,j)+Vp(i-1,j))+((dx/dy)^2).*((Vp(i,j+1)+Vp(i,j-1))))/(2*(1+(dx/dy)^2));
end;


%%
%Plot output from numerical solution
figure(1);
imagesc(x,y,V');
axis xy;
axis image;
colorbar;
colormap(gray(1024));
xlabel('X'); ylabel('Y');
title('Numerical Solution to Laplace Equation');

   