% Program 8.5 Finite difference solver for 2D Poisson equation 
%   with Dirichlet boundary conditions on a rectangle 
% Input: rectangle domain [xl,xr]x[yb,yt] with MxN space steps 
% Output: matrix w holding solution values 
% Example usage: w=poisson(0,1,1,2,4,4) 
function w=poisson7(xl,xr,yb,yt,M,N)  
K=1.68;
delta=0.1;
P=5;
L=2;
m=M+1;n=N+1; mn=m*n; 
h=(xr-xl)/M;h2=h^2;k=(yt-yb)/N;k2=k^2; 
x=xl+(0:M)*h; % set mesh values 
y=yb+(0:N)*k; 
A=zeros(mn,mn);b=zeros(mn,1); 
H=0.005;
for i=2:m-1 % interior points 
  for j=2:n-1 
    A(i+(j-1)*m,i-1+(j-1)*m)=1/h2;
    A(i+(j-1)*m,i+1+(j-1)*m)=1/h2; 
    A(i+(j-1)*m,i+(j-1)*m)=-2/h2-2/k2-(2*H)/(K*delta); 
    A(i+(j-1)*m,i+(j-2)*m)=1/k2;
    A(i+(j-1)*m,i+j*m)=1/k2; 
    b(i+(j-1)*m)=0; 
  end 
end 
for i=2:m-1 % bottom and top boundary points    
  A(i,i+(0)*m)=-3/(2*k)+H/K;
  A(i,i+(1)*m)=4/(2*k);
  A(i,i+(2)*m)=-1/(2*k);
  b(i)=0; 
  A(i+(n-1)*m,i+(n-1)*m)=-3/(2*k)+H/K;
  A(i+(n-1)*m,i+(n-2)*m)=4/(2*k);
  A(i+(n-1)*m,i+(n-3)*m)=-1/(2*k);
  b(i+(n-1)*m)=0; 
end 
for j=1:n % left and right boundary points 
  A(1+(j-1)*m,1+(j-1)*m)=-3/(2*h)+H/K;
  A(1+(j-1)*m,2+(j-1)*m)=4/(2*h);
  A(1+(j-1)*m,3+(j-1)*m)=-1/(2*h);
  if j<=n/2
    b(1+(j-1)*m)=-P/(L*delta*K); 
  else
    b(1+(j-1)*m)=0; 
  end
  A(m+(j-1)*m,m+(j-1)*m)=-3/(2*h)+H/K;
  A(m+(j-1)*m,(m-1)+(j-1)*m)=4/(2*h);
  A(m+(j-1)*m,(m-2)+(j-1)*m)=-1/(2*h);
  b(m+(j-1)*m)=0; 
  if j>=n/2
    b(m+(j-1)*m)=-P/(L*delta*K);
  else
    b(m+(j-1)*m)=0;
  end
end 
v=A\b; % solve for solution in v labeling 
v=v+20;
w=reshape(v(1:mn),m,n); %translate from v to w
surf(x,y,w') 