% 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=cool(0,1,1,2,4,4) 
function w=cool(xl,xr,yb,yt,M,N) 
ub = 20; %degrees Celcius
K=1.68; %watts/cm degrees Celcius
H=.005; %watts/cm^2 degrees Celcius
delta=.1; %cm

P=5; %watts
L=2; %cm

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*delta); 
    A(i+(j-1)*m,i+(j-2)*m)=1/k2;         A(i+(j-1)*m,i+j*m)=1/k2; 
    
  end 
end 
for i=1:m % bottom and top boundary points 
  j=1;     A(i+(1-1)*m,i+(j-1)*m) = -3/(2*k) + H/K;      %bottom
  j=2;     A(i+(1-1)*m,i+(j-1)*m) = 4/(2*k);
  j=3;     A(i+(1-1)*m,i+(j-1)*m) = -1/(2*k);
  
  j=n-2;   A(i+(n-1)*m,i+(j-1)*m) =  -1/(2*k);      %top
  j=n-1;   A(i+(n-1)*m,i+(j-1)*m) = 4/(2*k);
  j=n;     A(i+(n-1)*m,i+(j-1)*m) =  -3/(2*k) + H/K;
  
  
end 
for j=2:n-1 % left and right boundary points 
  i=1;     A(1+(j-1)*m,i+(j-1)*m)=-3/(2*h);      b(i+(j-1)*m)= -P/(L*delta*K); %left
  i=2;     A(1+(j-1)*m,i+(j-1)*m)=4/(2*h);      %b(i+(j-1)*m)= -P/(L*delta*K);
  i=3;     A(1+(j-1)*m,i+(j-1)*m)=-1/(2*h);      %b(i+(j-1)*m)= -P/(L*delta*K);
  
  i=m-2;   A(m+(j-1)*m,i+(j-1)*m)= -1/(2*h);     %right
  i=m-1;   A(m+(j-1)*m,i+(j-1)*m)=4/(2*h);       
  i=m;     A(m+(j-1)*m,i+(j-1)*m)=  -3/(2*h) + H/K;    
end 
v=A\b; % solve for solution in v labeling 
v = v+ ub;
w=reshape(v(1:mn),m,n); %translate from v to w 
mesh(x,y,w')