%Problem 5.36 - fin cooling problem with constant x-sectional area
  %version with one-sided deivative b.c.
  clear;
  ts=100;tf=30;
  nn=7;n=nn-2;h=35;k=386;l=.03;dx=l/(nn-1);
  x=0:dx:l;xi=dx/2:dx:l;
  dia=0.0025;area=pi*dia^2/4*ones(nn,1);
  p=pi*dia*ones(nn,1);
  m2=h*p./(k*area);
  %set far end boundary condition febc => 0 for adiabatic or 1 for convective
  febc=1;
  %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*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)
  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(' Total number of fins')
  nfins=1*1/(.006*.006)
  disp(' Total Heat Flow rate for whole Fin in Watts')
  totq=q*nfins+h*(1-nfins*area(1))*(ts-tf)
  disp(' Fin Efficiency (%) ')
  %calculate fin efficiency
  surfarea=trapz(x,p);
  effn=q/(h*surfarea*t(1))*100
  disp(' Fin Effectiveness (%) ')
  %calculate fin effectiveness
  effc=q/(h*area(1)*t(1))*100
  % plot theta vs x
  plot(x,temp)