% Reality Check 8 Heat distribution on a cooling fin 
%   with Robinson 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=fin2(0,4,0,4,10,10)
% Example usage: w=fin2(0,4,0,4,50,50)
% Example usage: w=fin2(0,4,0,4,100,100)

function w=fin2(xl,xr,yb,yt,M,N)
format long
f=@(x,y) 0; % define input function data 

% Constants
H=0.005; K=1.68; ub=20; d=0.1; L=2; P=5;

% Intervals 
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; 

% matrix A & b 
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)=f(x(i),y(j)); 
  end 
end 
for i=1:m % bottom and top boundary points 
  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; 
  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/2 % power enters along the left side
  i=1; 
  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); 
end
for j=(N/2)+1:n-1 % left boundary points 
  i=1; 
  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
for j=2:n-1 % right boundary points 
  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+ub; % solve for solution in v labeling 
w=reshape(v(1:mn),m,n); %translate from v to w 
mesh(x,y,w')
view(35,35);
xlabel('x'); ylabel('y'); zlabel('Temp.');
maxT=max(v)