function y = f(x)

%---------------
% Ca    -> x(1)
% CO3   -> x(2)
% H     -> x(3)
% HCO3  -> x(4)
% H2CO3 -> x(5)
%---------------

K1 = 10^-7.1;
K2 = 10^-10.25;
K3 = 10^-6.37;
K4 = 5.96e-9;
Kw = 10^-14;
OH = 10^-8;
Pco2 = 10^-3.5;

y = zeros(size(x));

y(1) = K1*Pco2 - x(3)*x(4);
y(2) = K2*x(4) - x(3)*x(2);
y(3) = K3*x(5) - x(4)*x(3);
y(4) = K4 - x(1)*x(2);
y(5) = x(3) - x(4) - x(2) + x(5)+OH;
y(6) = x(1) - 1e-3;
y(7) = Kw-x(3)*OH;

