% Separation of variables solution to transient conduction in simple geometries (Table 6-1) 
function [T, Bi, tt, xx, TT] = TSepVar(domain,x,t,varargin)
  [id,a,k,L,h,Tinf,To]=Unpack(domain);
  if length(varargin) N=varargin{1};
  else N=50;
  end
  Bi=h*L/k;
  tt=a*t(:)'/L^2;
  xx=x/L;
  if id=='Rec'
    EigenErr = @(b)b.*sin(b)-Bi*cos(b);
    Bn=EigenVal(EigenErr,N,pi);
    Cn = 4*sin(Bn)./(2*Bn+sin(2*Bn));
    EigenFun = cos(Bn*xx(:)');
  elseif id=='Cyl'
    EigenErr = @(b)b.*besselj(1,b)-Bi*besselj(0,b);
    Bn=EigenVal(EigenErr,N,pi);
    Cn = 2*besselj(1,Bn)./Bn./(besselj(0,Bn).^2+besselj(1,Bn).^2);
    EigenFun = besselj(0,Bn*xx(:)');
  elseif id=='Sph'
    xx(find(xx==0))=1e-6; % no zeros
    EigenErr = @(b)(1-Bi)*sin(b)-b.*cos(b);
    Bn=EigenVal(EigenErr,N,pi,1); % flag to exclude 0 from roots 
    Cn = 4*(sin(Bn)-Bn.*cos(Bn))./(2*Bn-sin(2*Bn));
    Br=Bn*xx(:)';
    EigenFun = sin(Br)./(Br);
  else
    error('domain id must = ''Rec'', ''Cyl'' or ''Sph''');
  end
  Bn2 = Bn.*Bn;
  TT = ( Cn.*exp(-Bn2*tt) )' * EigenFun;
  T=Tinf+TT*(To-Tinf);
end
