function wd()

global MO

PhysicalConstants();

rho_c = 10.^(3:0.3:12); % (g/cm^3) array of central densities from which
            % we start integrations for the mass versus radius p[lot

%rho_c = 10.^(3:1:5); % (g/cm^3) array of central densities from which

nr = 1000; % number of points in the span of radii for DE solver
rmax = 1e10; % (cm) max radius to which we attempt to integrate
rmin = 1e5; % (cm) min radius from which we start the integrartion

% Plot radius versus mass to see Chandrasekhar limit
for ir = 1:length(rho_c)
[r{ir}, m{ir}, rho{ir}, rx(ir), mx(ir)] = OneProfle(rho_c(ir)...
   , nr, rmin, rmax);
end
f101 = figure(101);
% mass in solar masses, radius in 1000 km
plot( mx/MO, rx/1e8, ':.', 'LineWidth', 2);
ylabel('WD radius/(1000 km)','FontSize',12, 'FontWeight','bold');
xlabel('WD mass/Msun','FontSize',12, 'FontWeight','bold');
set(gca, 'FontSize',12, 'FontWeight','bold');
print(f101,'MassRadius','-dpng')

% Plot interior properties for special cases of 0.6 and 1.2 solar masses
% To get this we start with log10(rho_c) = 6.5 and 8.1 respectively
[r06, m06, rho06, rx06, mx06] = OneProfle(10^6.5, nr, rmin, rmax);
p06 = Pressure(rho06);
[r12, m12, rho12, rx12, mx12] = OneProfle(10^8.1, nr, rmin, rmax);
p12 = Pressure(rho12);
f102 = figure(102);
subplot(3,1,1)
plot( r06/1e8, m06/MO, r12/1e8, m12/MO, 'LineWidth', 2);
ylabel('m(r) / M Sun','FontSize',12, 'FontWeight','bold');
legend('M/Sun=0.6', 'M/Sun=1.2')
subplot(3,1,2)
plot( r06/1e8, log10(p06), r12/1e8, log10(p12), 'LineWidth', 2);
ylabel('log10(P(Pa))','FontSize',12, 'FontWeight','bold');
subplot(3,1,3)
plot( r06/1e8, log10(rho06), r12/1e8, log10(rho12), 'LineWidth', 2);
ylabel('log10(\rho(g/cm^3))','FontSize',12, 'FontWeight','bold');
xlabel('r/(1000 km)','FontSize',12, 'FontWeight','bold');
print(f102,'InternProp','-dpng')

% Effective polytropic index
P = Pressure(rho_c);
dlnPdlnrho = rho_c./P.*pressureDerivative(rho_c);

Problem 6. Nuclear fusion rate A hot plasma (kT = 10 keV) consis...