%Problem - fin cooling problem with variable x-sectional area
  %version with one-sided deivative b.c.
  %time-dependent version - implicit method
  clear;
  nn=11;n=nn-2;h=20;k=20;l=2/12;dx=l/(nn-1);
  x=0:dx:l;xi=dx/2:dx:l;
  area=.5*(1-1/3*sinh(x/l))/144;p=sqrt(area);
  m2=h*p./(k*area);
  tbase=120;
  dt=0.01;alf=1.04;lbd=alf*dt/dx^2;% ***** added or modified for time solution
  %set far end boundary condition febc => 0 for adiabatic or 1 for convective
  febc=0;
  %tri-diagonal coeffients
  for i=1:n
     a(i)=lbd*(area(i)+area(i+1))/(2*area(i+1));% ***** added or modified for time solution
     c(i)=lbd*(area(i+1)+area(i+2))/(2*area(i+1));% ***** added or modified for time solution
  end
  b=-1+lbd*(-2-m2(2:nn-1)*dx*dx);% ***** added or modified for time solution
  %initial conditions
  t=tbase*diag(eye(nn));d=-t(2:nn-1);% ***** added or modified for time solution 
  %modify coeffts for 2nd-order b.c. - elimination
  d(1)=d(1)-a(1)*t(1); denom=(3/2+febc*h*dx/k);% ***** added or modified for time solution
  a(n)=a(n)-.5*c(n)/denom;b(n)=b(n)+2*c(n)/denom;
  %time dependent part
  for j=1:20% ***** added or modified for time solution
     time=j*dt
     %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 new rhs
     d(2:nn-1)=-t(3:nn);d(1)=-a(1)*t(1)-t(2);% ***** added or modified for time solution
  end
  disp(' Heat Flow rate at Base of Fin')
  %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,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,t)
