%
%   Fuel_Temp_SH1.M   Solve for T(r) in a cylindrical fuel pin (kf = const)
%                        Analytically and via the Shooting Method
%
%   Given differential eqn.:   T'' + (1/r)T' + qqq/kf =  0   
%                        with  T'(0) = 0  and   T(Rf) = Ts 
%
%
%   File prepared by J. R. White, UMass-Lowell (Sept. 2003)
%

%
%   getting started
      clear all, close all,  nfig = 0;
%
%   set properties
      global qqq kf
      qqq = 1.35e+7;     % average volumetric heat source (BTU/hr-ft^3)
      kf = 1.1;          % thermal conductivity of UO2 (BTU/hr-ft-F)
      Rf = 0.0175;       % fuel pin radius (ft)
      Ts = 650;          % surface temperature (F)
%
%   Analytical Solution (exact)
%
      re = linspace(0,Rf,30);
      Te = Ts + (qqq/(4*kf))*(Rf*Rf - re.*re);
%
%   Solution via Shooting Method (numerical approx)
%
%   set coefficients for BCs and domain for independent variable
      xo = 0;  xf = Rf;
      zbc = [0 1.0 0.0;      % left BC  --> T'(xo) = 0.0
             1 0 Ts];       % right BC --> T(xf) = Ts
%
%   solution via Shooting Method (numerical approx)
      tol = 1e-6;  options = odeset('RelTol',tol);   % set tight tolerance for ODE soln
      [xs,zs] = bvp2sh('fuel_temp_sh1_eqns',[xo xf],zbc,options);
%
%   plot solution 
      nfig = nfig+1;   figure(nfig)
      plot(re,Te,'r-',xs,zs(:,1),'bo','LineWidth',2),grid
      title('Fuel Pin Radial Temperature Profile')                  
      xlabel('Radial Distance (ft)'),ylabel('Temperature (F)')    
      legend('Analytical Solution','Shooting Method')                    
%
%   end program
%
