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, indexAtCrossectionY] = min(abs(logcNP-crossectionY));
crossectionX = pH(indexAtCrossectionY(1))
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
%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;
.
