Obtained weird answers to my Numerical Integration

I've used I_s = 82 and C = 0.009 but my answers obtained are in the range of millions. Any help on this coding?

[SCd indent code]

Vol_liq = 1.5*10^-3;
Q = 6.5*10^-6;
d_r = 0.08;
d = 0.1;

%Calculate and display cross-sectional area of reactor
A = (pi/4)*d^2;

%Calculate and display cross-sectional area of riser
A_r = (pi/4)*d_r^2;

%Calculate and display cross-sectional area of downcomer
A_d = A-A_r;

%Calculate and display superficial gas velocity
U = Q/A_r;

%Calculate and display Liquid circulation velocity [Bello et al.(1985)]
u_lr = 0.66*U^(0.33)*(A_d/A_r)^(0.78);

%Prompt user for apparent viscosity to calculate and display gas holdup [Popovic and Robinson,1984] in reactor
n = 10^-3;
e_g = 0.465*U^(0.65)*(1+A_d/A_r)^(-1.06)*n^(-0.103);

%Calculate and display linear liquid velocity in riser
J_lr = u_lr/(1-e_g);

%Calculate and display liquid flow rate in reactor
Q_l = J_lr*A_r*(1-e_g);

%Calculate and display linear liquid velocity in downcomer
J_ld = Q_l/A_d;

%Prompt user for height of riser(draft tube), reactor, bottom
H = 0.3;
H_d = 0.11;
H_b = 0.04;

%Calculate and display height of separator section
H_s = H*(1+e_g)-H_d-H_b;

%Calculate and display residence time for riser,downcomer and separator
%regions of bioreactor
t_r = (H_d+H_b)/J_lr;
t_d = (H_d+H_b)/J_ld;
t_s = (H_s*A*(1-e_g))/Q_l;
t_cycle = t_r + t_d + t_s;

Time_riser = ['Residence time in riser section = ' num2str(t_r) ' s'];
Time_downcomer = ['Residence time in downcomer section = ' num2str(t_d) ' s'];
Time_separator = ['Residence time in separator section = ' num2str(t_s) ' s'];
Time_cycle = ['Total time for 1 cycle of circulation = ' num2str(t_cycle) ' s']; 

disp(Time_riser);
disp(Time_downcomer);
disp(Time_separator);
disp(Time_cycle);

%-------------------------------------------------------------------------------------------------------------------------
%Prompts the user to input value of the intensity of light source and initial concentration of cell
I_s = 82;
Cell = 0.009*1.3*10^8;

r_riser = linspace(0,40,40);
r_downcomer = linspace(42,50,8);
r_separator = linspace(0,50,50);

I_riser = 2.*pi.*r_riser.*I_s.*(exp(-0.079.*(50-r_riser))+6*10^(-9)*Cell^(-2.5))./(exp(3.5*10^(-10).*Cell.*(50-r_riser))+6*10^(-9)*Cell^(-2.5));
I_ave_riser = trapz(I_riser,r_riser)/A_r;

I_downcomer = 2.*pi.*r_downcomer.*I_s.*(exp(-0.079.*(50-r_downcomer))+6*10^(-9)*Cell^(-2.5))./(exp(3.5*10^(-10).*Cell.*(50-r_downcomer))+6*10^(-9)*Cell^(-2.5));
I_ave_downcomer = trapz(I_downcomer,r_downcomer)/A_d;

I_separator = 2.*pi.*r_separator.*I_s.*(exp(-0.079.*(50-r_separator))+6*10^(-9)*Cell^(-2.5))./(exp(3.5*10^(-10).*Cell.*(50-r_separator))+6*10^(-9)*Cell^(-2.5));
I_ave_separator = trapz(I_separator,r_separator)/A;

%Display average intensity of each section
fprintf('Average intensity in the riser = (%6.7f) umol/m^2-s\n',I_ave_riser);
fprintf('Average intensity in the downcomer = (%6.7f) umol/m^2-s\n',I_ave_downcomer);
fprintf('Average intensity in the separator section = (%6.7f) umol/m^2-s\n',I_ave_separator);

Answer 1

[Sean de Wolski]

Use the debugger (set a break point on the first line) and step through line by line. Inspect the value of the variables at each line. When do they explode? You'll be able to answer your own question and learn a fundamental of programming all at once!

[Lu]

It exploded during the numerical integration and I have checked the previous lines though....

[Lu]

I've tried to debug the trapezoidal integration but it gave me an error: Subscript indices must either be real positive integers or logicals.

[Jan]

Then I suggest to fix this problem.