%Problem 5.25 - fin cooling problem with variable x-sectional area
  %version with one-sided derivative b.c.
  clear;
  ts=180;tf=25;
  nn=21;n=nn-2;hc=25;k=180;l=.05;dx=l/(nn-1);
  emiss=0.9;stefanboltz=5.6697e-8;
  x=0:dx:l;xi=dx/2:dx:l;
  wd=1;bd=.01*(1-x/l);area=wd*bd;
  p=2*(wd+bd);
  %initially set radiation heat transfer to zero - recalculate after
  %temperature is calculated below.
  hr=zeros(1,nn);
  %Overall heat transfer coefficient will be non-linear if h due to radiation
  %added. Therefore iteration will be necessary.
  for iter=1:4
   h=hr+hc;   
  m2=h.*p./(k*area);
  %set far end boundary condition febc => 0 for adiabatic or 1 for convective
  febc=0;
  %tri-diagonal coeffients
  for i=1:n
     a(i)=(area(i)+area(i+1))/(2*area(i+1));
     c(i)=(area(i+1)+area(i+2))/(2*area(i+1));
     end
  b=(-2-m2(2:nn-1)*dx*dx);
  d=diag(zeros(n)); 
  %modify coeffts for 2nd-order b.c. - elimination
  t(1)=ts-tf; 
  d(1)=-a(1)*t(1); denom=(3/2+febc*h(n)*dx/k);
  a(n)=a(n)-.5*c(n)/denom;b(n)=b(n)+2*c(n)/denom;
  %call tri-diagonal solver
  disp(' Temperature (theta) above ambient')
  t(2:nn-1)=trid(a,b,c,d);
  t(nn)=-.5/denom*t(nn-2)+2/denom*t(nn-1)
  %Calculate radiation heat transfer coefficient
  hr=emiss*stefanboltz*((t+tf+273.15).^2+290^2).*(t+tf+273.15+290);
  h=hr+hc;
  end
  disp(' Temperature T in degrees Centrigade')
  temp=t+tf
  disp(' Temperature at the tip in degrees Centrigade')
  temp(nn)
  disp(' Heat Flow rate at Base of Fin in Watts')
  %calculate heat flow rate at x=0
  dtdx=(1.5*t(3)-2*t(2)+.5*t(1))/dx;
  q=-k*area(1)*dtdx
  disp(' Fin Efficiency (%) ')
  %calculate fin efficiency
  surfarea=trapz(x,h.*p);
  effn=q/(surfarea*t(1))*100
  disp(' Fin Effectiveness (%) ')
  %calculate fin effectiveness
  effc=q/(h(1)*area(1)*t(1))*100
  % plot theta vs x
  plot(x,temp)
