MATLAB takes too long time to find solution

Hello everyone, I'm ChemE student trying to solve what is called %recovery of two gases in a stripper column. there's bisection method inside the while loop and the while loop has two conditions. the problem solved when the value of 'er_h2s' and 'er_co2' are lower than 'tol'

When I ran the m-file, there was no sign of working but there was a pause button on the top tab. I was certain it was doing its job, I was waiting for about 20 minutes and still nothing popped up on the screen. here's my codes: (I also attach them)

clc;

%nomor tray sesungguhnya = i-1; contoh tray 1 = (2-1);
q = 252.9745;
lambda = 2202.59;
P = 2*101325;
i = 2; %COBA 2 dulu. angka 2 artinya tray 1
tol = 1e-5;
K1 = 0.3354;
K2 = 0.5766;
K4 = 18.0273;
K5 = 1.97433e-07;
K6 = 5.091610747;
K7 = 3.3986e-12;

n_MDEAH(1) = 0.02158;
n_PZH(1) = 1.15195e-9;
n_HCO3(1) = 0.020554338;
n_HS(1) = 0.001027717;
n_MDEA(1) = 0.026672116;
n_W(1) = 0.348588637;
n_PZ(1) = 0.004450412;
n_CO3(1) = 3.09926e-5;
n_H(1) = 3.84719e-8;
n_OH(1) = 5.92995e-11;
n_CH4(1) = 0.8316;

He_H2S = 4.71e12;
He_CO2 = 1.54e15;

er_h2s(2) = 1;
er_co2(2) = 1;

while er_h2s(i) > tol && er_co2(i) > tol
rec_H2S_tray(i) = 0.4;
rec_CO2_tray(i) = 0.09;
n_H2S_out(i) = rec_H2S_tray(i)*n_HS(i-1);
n_CO2_out(i) = rec_CO2_tray(i)*n_HCO3(i-1);
n_HS_reacted(i) = n_H2S_out(i-1);
n_HCO3_reacted(i) = n_CO2_out(i-1);
    
%~~~~START EQUILIBRIUM CALCULATION USING BISECTION~~~~
x1=0.0044; %guess awal
n_PZH(i) = (n_PZ(i-1) + n_PZH(i-1)) - x1;
n_MDEAH(i) = (n_MDEAH(i-1) - (n_HCO3_reacted(i) - n_PZH(i-1)) - n_HS_reacted(i)) - n_PZH(i);
n_MDEA(i) = (n_MDEAH(i-1) + (n_HCO3_reacted(i) - n_PZH(i-1)) + n_HS_reacted(i)) + n_PZH(i);
n_OH(i) = K6*n_MDEA(i)/n_MDEAH(i);
n_H(i) = n_PZH(i)/(K6*x1);
f1 = K7*n_H(i)*n_OH(i);
 
x2=0.01; %guess awal
n_PZH(i) = (n_PZ(i-1) + n_PZH(i-1)) - x2;
n_MDEAH(i) = (n_MDEAH(i-1) - (n_HCO3_reacted(i) - n_PZH(i-1)) - n_HS_reacted(i)) - n_PZH(i);
n_MDEA(i) = (n_MDEAH(i-1) + (n_HCO3_reacted(i) - n_PZH(i-1)) + n_HS_reacted(i)) + n_PZH(i);
n_OH(i) = K6*n_MDEA(i)/n_MDEAH(i);
n_H(i) = n_PZH(i)/(K6*x2);
f2=K7*n_H(i)*n_OH(i);

while f1*f2>=0
x1=input('masukkan nilai x1 = ' );
x2=input('masukkan nilai x2 = ' );
n_PZH(i) = (n_PZ(i-1) + n_PZH(i-1)) - x1;
n_MDEAH(i) = (n_MDEAH(i-1) - (n_HCO3_reacted(i) - n_PZH(i-1)) - n_HS_reacted(i)) - n_PZH(i);
n_MDEA(i) = (n_MDEAH(i-1) + (n_HCO3_reacted(i) - n_PZH(i-1)) + n_HS_reacted(i)) + n_PZH(i);
n_OH(i) = K6*n_MDEA(i)/n_MDEAH(i);
n_H(i) = n_PZH(i)/(K6*x1);
f1=K7*n_H(i)*n_OH(i);
    
n_PZH(i) = (n_PZ(i-1) + n_PZH(i-1)) - x2;
n_MDEAH(i) = (n_MDEAH(i-1) - (n_HCO3_reacted(i) - n_PZH(i-1)) - n_HS_reacted(i)) - n_PZH(i);
n_MDEA(i) = (n_MDEAH(i-1) + (n_HCO3_reacted(i) - n_PZH(i-1)) + n_HS_reacted(i)) + n_PZH(i);
n_OH(i) = K6*n_MDEA(i)/n_MDEAH(i);
n_H(i) = n_PZH(i)/(K6*x2);
f1=K7*n_H(i)*n_OH(i);
end

e=1;
ite=0;

while e>=tol
x3=(x1+x2)/2;
n_PZH(i) = (n_PZ(i-1) + n_PZH(i-1)) - x3;
n_MDEAH(i) = (n_MDEAH(i-1) - (n_HCO3_reacted(i) - n_PZH(i-1)) - n_HS_reacted(i)) - n_PZH(i);
n_MDEA(i) = (n_MDEAH(i-1) + (n_HCO3_reacted(i) - n_PZH(i-1)) + n_HS_reacted(i)) + n_PZH(i);
n_OH(i) = K6*n_MDEA(i)/n_MDEAH(i);
n_H(i) = n_PZH(i)/(K6*x3);
f3=K7*n_H(i)*n_OH(i);    

ite=ite+1;
e=abs(x1-x2)/2;

if f1*f3<=0
    x2=x3;
    f2=f3;
else
    x1=x3;
    f1=f3;
end

end
%~~~~END OF EQUILIBRIUM CALCULATION USING BISECTION~~~~

 n_PZ(i)= x3;
 n_HS(i) = n_HS(i-1) - n_HS_reacted(i);
 n_HCO3(i) = n_HCO3(i-1) - n_HCO3_reacted(i);
 n_CO3(i) = K4*n_HCO3(i)*n_OH(i);
 n_H2S(i) = n_MDEAH(i)*n_HS(i)/(K2*n_MDEA(i));
 n_CO2(i) = n_MDEAH(i)*n_HCO3(i)/(K1*n_MDEA(i));
 
 G(i+1) = n_H2S_out(i)+n_CO2_out(i)+n_CH4(1)+0.85*q/lambda;
 y_H2S(i) = n_H2S_out(i)/G(i+1);
 y_CO2(i) = n_CO2_out(i)/G(i+1);
 
 %~~~INTERFACIAL CONCENTRATION CALCULATION~~~
 %[A]i = yA.P/HeA
 Mi_H2S(i) = y_H2S(i)*P/He_H2S;
 Mi_CO2(i) = y_CO2(i)*P/He_CO2;
 
 %Difusivitas
 D_CO2MDEA = 9.675e-8;
 D_MDEAW = 1.2046e-7;
 D_H2SW = 9.5384e-7;
 
 %Mass transfer coefficient
 kL_CO2 = 5.381e-5;
 kL_H2S = 6.2624e-5;
 
 %~~~ENHANCEMENT FACTOR~~~
 k_CO2 = 3.21e-9;
 k_PZ = 272739.3;
 k1 = k_CO2*n_MDEA(i)+k_PZ*n_PZ(i);
 Hatta_CO2 = D_CO2MDEA*n_MDEA(i)*k1/(kL_CO2^2);
 E_CO2 = (Hatta_CO2)^0.5;
 
 k_H2S = 5.52e-9;
 k2 = k_H2S*n_MDEA(i);
 E_H2S = 1 + n_MDEA(1)*D_MDEAW/(Mi_H2S(i)*D_H2SW);
 
 %~~~~~~MASS TRANSFER FLUX CALCULATION~~~~~~
 N_H2S(i) = kL_H2S*E_H2S*(Mi_H2S(i) - n_H2S(i));
 N_CO2(i) = kL_CO2*E_CO2*(Mi_CO2(i) - n_CO2(i));
 
 %~~~~~Liquid mass balance~~~
 a = 0.007184;
 n_HS_flux(i) = n_HS(1) - N_H2S(i)*a;
 n_HCO3_flux(i) = n_HCO3(1) - N_CO2(i)*a;
 
 %~~~~~new recovery calculation~~~~
 rec_H2S_tray_new(i) = (n_HS(i-1) - n_HS_flux(i))/n_HS_flux(i);
 rec_CO2_tray_new(i) = (n_HCO3(i-1) - n_HCO3_flux(i))/n_HCO3_flux(i);
 
 er_H2S = abs((rec_H2S_tray(i) - rec_H2S_tray_new(i))/rec_H2S_tray_new(i));
 er_CO2 = abs((rec_CO2_tray(i) - rec_CO2_tray_new(i))/rec_CO2_tray_new(i));
 
 rec_H2S_tray(i) = rec_H2S_tray_new(i);
 rec_CO2_tray(i) = rec_CO2_tray_new(i);
end

disp(['%recovery H2S tray =',num2str(rec_H2S_tray(i))]);
disp(['%recovery H2S tray new =',num2str(rec_H2S_tray_new(i))]);
disp(['ite= ', num2str(ite)]);

I have no idea of what's going on, maybe there's something wrong with my codes? Please help.. Thank you

[KALYAN ACHARJYA]

Here multiple while loops, till the condition is true, the loop goes within.

Answer 1

[per isakson]

See Profile to Improve Performance

I let your code run for a minute with the profiler active.

There are a lot of warnings like

The variable 'n_H2S_out' appears to change size on every loop iteration. 
Consider preallocating for speed.

It converts very slowly

[Bertiningrum Devi]

Thank you for your response, Per.

my goodness..it looks terrible

how do I overcome those problems? I need some help..

[per isakson]

There a couple of minor problems

• scalars are refered to with an index, e.g. n_CH4(1) = 0.8316; If nothing else it's confusing. Replace n_CH4(1) by n_CH4, et cetera.
• it's vectors with the length 2, e.g. rec_H2S_tray, which causes all the messages, Consider preallocating for speed. Add rec_H2S_tray = nan(1,2); at the top of the code. It will save a bit of time and make the messages disappear.

Then there is the large problem; the code converges slowly if at all. I don't understand what numerical method hides in your code. I haven't really tried. First question:

• are you sure that the code does what you believe it to do?

I added a couple of print statements

er_co2(2) = 1;
while er_h2s(i) > tol && er_co2(i) > tol
    fprintf( '\ner_h2s = %12f, %12f\n', er_h2s(i), er_co2(i) );  % <<<<<<
    rec_H2S_tray(i) = 0.4;

and

f2=K7*n_H(i)*n_OH(i);
disp(' ')    % <<<<<<<<<<<<<<
while abs(f1*f2)>=0
    fprintf( 'f1,f2 = %12f, %12f\n', f1, f2 ); % <<<<<<<<<<<<<<
    x1=input('masukkan nilai x1 = ' );         % ??????????????
    x2=input('masukkan nilai x2 = ' );         % ??????????????    

I added abs() because the product was negative. Execution never entered into this while-loop. The input()-statement makes the code impractical to run. You need to assign guesses automatically.

And a last print-statement.

ite=0;
disp(' ')  % <<<<<<<<<<<<<<
while e>=tol
    fprintf( 'e = %f\n', e ); % <<<<<<<<<
    x3=(x1+x2)/2;

Start the code and click Pause to halt the execution. The command window is full of output, which you should try to interpret.

And finally see: Debug a MATLAB Program and Examine Values While Debugging