% Clear variables and figures
clear;
clf;
global h;
H=380;
% Problem 1.d find the position h from the criticality eqn.
h = fzero(@criticalityCondition,175);
disp(['Control rod position= ',num2str(h),' cm']);
% 2.b compute the flux profiles
zz1 = 0:1:h;
zz2 = h:1:H;
phi1 = flux1(zz1);
phi2 = flux2(h,zz2);
% plot the results
zz = [zz1,zz2];
phi = [phi1,phi2];
plot(zz,phi);
% Problem 2.c integrate the functions
first_integral = quad(@flux1,0,h);
second_integral = quad(@flux2s,h,H);
total_integral = first_integral + second_integral;
disp(['1st integral= ',num2str(first_integral)]);
disp(['2nd integral= ',num2str(second_integral)]);
disp(['Total integral= ',num2str(total_integral)]);



