Untitled

Kh=10^(-1.46);
K1=10^(-6.3);
K2=10^(-10.3);
Kw=10^(-14);
Pco2=@(t) 0.011825.*(t-1980.5).^2+1.356975.*(t-1980.5)+339;
t=[1958 2003];
Pco2_1958= Pco2(t(1));
Pco2_2003= Pco2(t(2));
disp(['Pco2 1958:', num2str(Pco2_1958)]);
disp(['Pco2 2003:', num2str(Pco2_2003)]);

fun_1958= @(H) K1./(1e6.*H)*Kh*Pco2_1958+2*K2*K1./(1e6.*H.^2)*Kh*Pco2_1958+Kw./H-H;
fun_2003= @(H) K1./(1e6.*H)*Kh*Pco2_2003+2*K2*K1./(1e6.*H.^2)*Kh*Pco2_2003+Kw./H-H;
Kh=10^(-1.46);
K1=10^(-6.3);
K2=10^(-10.3);
Kw=10^(-14);
H0=[10^(-7) 10^(-7)];
concentrazione_1958= fsolve(fun_1958,H0);
concentrazione_2003= fsolve(fun_2003,H0);
disp(['C 1958:', num2str(concentrazione_1958)]);
disp(['C 2003:', num2str(concentrazione_2003)]);

pH_1958= -log10(concentrazione_1958);
pH_2003= -log10(concentrazione_2003);
disp(['pH 1958:', num2str(pH_1958)]);
disp(['pH 2003:', num2str(pH_2003)]);