function w=poisson2(M,N,P)
xl = 0;
xr = 4;
pl = 0;
pr = 2;
yb = 0;
yt = 4;
H = 0.005; K = 1.68; d = 0.1; 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); 
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*d)); 
    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=1:m % bottom and top boundary points 
  %bottom
  j=1;
  A(i+(j-1)*m,i+(j-1)*m)= -3/(2*k) + H/K;
  A(i+(j-1)*m,i+j*m)=4/(2*k);
  A(i+(j-1)*m,i+(j+1)*m)=-1/(2*k);
  b(i+(j-1)*m)=0; 
  %top
  j=n;
  A(i+(j-1)*m,i+(j-1)*m)=-3/(2*k)+H/K;
  A(i+(j-1)*m,i+(j-2)*m)= 4/(2*k);    
  A(i+(j-1)*m,i+(j-3)*m)=-1/(2*k);
  b(i+(j-1)*m)=0; 
end 
for j=2:n-1 % left and right boundary points 
  %left
  i=1;
  if ((j/M)*(yt - yb) < pr && (j/M)*(yt - yb) >= pl)
      A(i+(j-1)*m,i+(j-1)*m)   =-3/(2*h);
      A(i+(j-1)*m,i+1+(j-1)*m) = 4/(2*h);
      A(i+(j-1)*m,i+2+(j-1)*m) = -1/(2*h);
      b(i+(j-1)*m)             =-P/((L)*d*K);
  else
      A(i+(j-1)*m,i+(j-1)*m) = -3/(2*h)+ H/K;
      A(i+(j-1)*m,i+1+(j-1)*m) = 4/(2*h);
      A(i+(j-1)*m,i+2+(j-1)*m) = -1/(2*h);
      b(i+(j-1)*m)=0;
  end
  %right
  i=m;
  A(i+(j-1)*m,i+(j-1)*m)=-3/(2*h)+H/K;
  A(i+(j-1)*m,i-1+(j-1)*m) = 4/(2*h);
  A(i+(j-1)*m,i-2+(j-1)*m) = -1/(2*h);
  b(i+(j-1)*m)=0; 
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') 
title(sprintf('Heat Distribution in Celsius on a 4 x 4 cm Fin with 1 mm Thickness, Grid Size %d x %d',M,N))
    xlabel('cm')
    ylabel('cm')
    zlabel('Temperature (C)') 
    view(65,15)
shading interp