Why my code find a second intersection of curve with line and not the first?

I am trying to find a intersection of very complex curve's equation with line parallel with X-axis, but Matlab gives me a result which is suitable for the second intersection. I need first one.

Plot code for better understand is also included. Magenta line should reach and stop at first intersection (pH=6.4 not the second one at pH=13,38). What am I doing wrong?

Oh yes, and I need a exact value of X-coordinate of the first intersection and also a graphical solution as it can be seen.

Here is a MWE:

clc;
close all;
%There are input constants
cP=0.05;
cN=0.05;
cMg=0.005;

logb0=2.88;
logb1=1.17;
logb2=4.85;
logb3=2.58;
logKs1=-12.6;

pKa1=-2.143; 
pKa2=-7.205;
pKa3=-12.34;

pKb1NH3=-4.753355133;

%there is a range of x-axis (0-14)
pH=0:0.01:14;


%There are calculations of cH2PO4 cHPO4 cPO4, they needs to be calculated for final curve equation (they are variables of final curve equation and depend on pH value)
Ka1=power(10,pKa1);
Ka2=power(10,pKa2);
Ka3=power(10,pKa3);

cH2PO4=cP.*Ka1.*(power(10,-pH)).^(2)./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
cHPO4=cP.*Ka1.*Ka2.*power(10,-pH)./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
cPO4=cP.*Ka1.*Ka2.*Ka3./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);

%There is calculations of cNH4, it have to be calculated for final curve equation (it is variable of final curve equation and depends on pH value also)
pOH=14-pH;
Kb1=power(10,pKb1NH3);
cNH4=cN.*Kb1./(power(10,-pOH)+Kb1);

%there are some partial calculations for the final curve equation
b0=power(10,logb0);
b1=power(10,logb1);
b2=power(10,logb2);
b3=power(10,logb3);
Ks1=power(10,logKs1);
cOH=power(10,pH-14);

%this is the final curve equation
logcNP=Ks1./(cPO4.*cNH4).*(1+b1.*cH2PO4+b0.*cHPO4+b2.*cPO4+b3*cOH);

%I am trying to calculate a intersection of the curve above with y value of 0.005 (cMg at the start of the code).
crossectionY=cMg

crossectionY = 0.0050

[crossectionY, indexAtCrossectionY] = min(abs(logcNP-crossectionY));
crossectionX = pH(indexAtCrossectionY(1))

crossectionX = 13.3800

%This is plot of the final curve equation (black) and calculated intersection with magenta line 
semilogy(pH,logcNP,'-', 'color', 'k', 'MarkerFaceColor','w', 'MarkerEdgeColor', 'm');
box off;
hold on;
plot([0 crossectionX],[cMg cMg],'-', 'color', 'm', 'LineWidth',0.5);
box off;
hold on;
axis([0 14 0.000000000000000001 1000000000000000000]);
xlabel('pH (-)');
ylabel('c {(mol/dm^{3})}');
axis square;
grid on;
axis on;
box off;
hold off;

Than you for answers. :-)

Answer 1

[Star Strider]

See if the lines oif copde I added do what you want —

%There are input constants
cP=0.05;
cN=0.05;
cMg=0.005;
logb0=2.88;
logb1=1.17;
logb2=4.85;
logb3=2.58;
logKs1=-12.6;
pKa1=-2.143; 
pKa2=-7.205;
pKa3=-12.34;
pKb1NH3=-4.753355133;
%there is a range of x-axis (0-14)
pH=0:0.01:14;
%There are calculations of cH2PO4 cHPO4 cPO4, they needs to be calculated for final curve equation (they are variables of final curve equation and depend on pH value)
Ka1=power(10,pKa1);
Ka2=power(10,pKa2);
Ka3=power(10,pKa3);
cH2PO4=cP.*Ka1.*(power(10,-pH)).^(2)./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
cHPO4=cP.*Ka1.*Ka2.*power(10,-pH)./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
cPO4=cP.*Ka1.*Ka2.*Ka3./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
%There is calculations of cNH4, it have to be calculated for final curve equation (it is variable of final curve equation and depends on pH value also)
pOH=14-pH;
Kb1=power(10,pKb1NH3);
cNH4=cN.*Kb1./(power(10,-pOH)+Kb1);
%there are some partial calculations for the final curve equation
b0=power(10,logb0);
b1=power(10,logb1);
b2=power(10,logb2);
b3=power(10,logb3);
Ks1=power(10,logKs1);
cOH=power(10,pH-14);
%this is the final curve equation
logcNP=Ks1./(cPO4.*cNH4).*(1+b1.*cH2PO4+b0.*cHPO4+b2.*cPO4+b3*cOH);
%I am trying to calculate a crosssection of the curve above with y value of 0.005 (cMg at the start of the code).
crossectionY=cMg

crossectionY = 0.0050

[crossectionY, indexAtCrossectionY] = min(abs(logcNP-crossectionY));
crossectionX = pH(indexAtCrossectionY(1))

crossectionX = 13.3800

xcidx = find(diff(sign(logcNP - 0.005)));                                   % ADDED: Approximate Crossing In dices For 'logcNP' At 0.005
for k = 1:numel(xcidx)
    idxrng = max(1,xcidx(k)-1) : min(numel(pH),xcidx(k)+1);                 % ADDED: Index Range For Interpolation
    pHx(k) = interp1(logcNP(idxrng),pH(idxrng), 0.005);                     % ADDED: Values Of 'pH' For 'logcNP = 0.005'
end

pHx

pHx = 1×2

6.3537 13.3819

%This is plot of the final curve equation (black) and calculated crossection with magenta line 
semilogy(pH,logcNP,'-', 'color', 'k', 'MarkerFaceColor','w', 'MarkerEdgeColor', 'm');
box off;
hold on;
plot([0 crossectionX],[cMg cMg],'-', 'color', 'm', 'LineWidth',0.5);
plot(pHx, ones(size(pHx))*0.005, 'sc', 'MarkerSize', 10)                    % ADDED: Plot Cyan Squares At Intersection Points
box off;
hold on;
axis([0 14 0.000000000000000001 1000000000000000000]);
xlabel('pH (-)');
ylabel('c {(mol/dm^{3})}');
axis square;
grid on;
axis on;
box off;
hold off;

.

[Josef]

Yeah thank you! It is doing what I need! :-)

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.