<html>
<head><title>coolingfin.m</title></head>
<body>
<pre>
% Program 8.5 Finite difference solver for Cooling fin equation
% with Robin boundary conditions on a rectangle
% Input: rectangle domain [xl,xr]x[yb,yt] with MxN space steps
%        P watts of power applied on left side from y=pb to pt
%        K conductivity of fin, H conductivity of surrounding medium
% Output: maximum heat across fin
% Example usage: coolingfin(0, 2, 0, 2, 0, 2, 5, 1.68, 0.005, 10, 10,50,150)
function w=coolingfin(xl,xr,yb,yt,pb,pt,P,K,H,M,N,plotmin,plotmax)
delta = 0.1; ub = 20;
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=sparse(mn,mn);b=zeros(mn,1);
idx=@(i,j) i + j * m - m; % expansion from 2D coords to 1D coord
for i=2:m-1 % interior points
    for j=2:n-1
        %                   uxx                +                   uyy                = 2*H/(K*width) * u
        % [u(i-1,j) - 2*u(i,j) + u(i+1,j)]/h^2 + [u(i,j-1) - 2*u(i,j) + u(i,j+1)]/k^2 - 2*H/(K*width) * u = 0
        A(idx(i,j), idx(i-1, j )) = 1/h2;
        A(idx(i,j), idx(i+1, j )) = 1/h2;
        A(idx(i,j), idx( i , j )) = -2/h2 - 2/k2 - 2*H/(K*delta);
        A(idx(i,j), idx( i ,j-1)) = 1/k2;
        A(idx(i,j), idx( i ,j+1)) = 1/k2;
        b(idx(i,j)) = 0;
    end
end
for i=1:m % bottom and top boundary points
    j=1;A(idx(i,j), idx(i,j))   = -H/K - 1.5/k;
	    A(idx(i,j), idx(i,j+1)) = 2/k;
		A(idx(i,j), idx(i,j+2)) = -0.5/k;
		b(idx(i,j)) = 0;
    j=n;A(idx(i,j), idx(i,j))   = -H/K - 1.5/k;
	    A(idx(i,j), idx(i,j-1)) = 2/k;
		A(idx(i,j), idx(i,j-2)) = -0.5/k;
		b(idx(i,j)) = 0;
end
for j=2:n-1 % left and right boundary points
    i=1;
	if yb+(j-1)*k >= pb && yb+(j-1)*k <= pt
	    A(idx(i,j), idx(i,  j)) = - 1.5/h;
	    A(idx(i,j), idx(i+1,j)) = 2/h;
		A(idx(i,j), idx(i+2,j)) = -0.5/h;
		b(idx(i,j)) = - P/((pt-pb)*delta*K);
	else
	    A(idx(i,j), idx(i,  j)) = -H/K - 1.5/h;
	    A(idx(i,j), idx(i+1,j)) = 2/h;
		A(idx(i,j), idx(i+2,j)) = -0.5/h;
		b(idx(i,j)) = 0;
	end
    i=m;A(idx(i,j), idx(i,  j)) = -H/K - 1.5/h;
	    A(idx(i,j), idx(i-1,j)) = 2/h;
		A(idx(i,j), idx(i-2,j)) = -0.5/h;
		b(idx(i,j)) = 0;
end
v=A\b; % solve for solution in v labeling
w=reshape(v(1:mn),m,n); %translate from v to w
mesh(x,y,w'+ub)
view([45, 45])
w=max(max(w))+ub;
    axis([xl xr yb yt plotmin plotmax])
    caxis([plotmin plotmax]);
end
</pre>
</body>