function y = f(x)

%----------------
% EDTAH3 -> x(1)
% OH     -> x(2)
% EDTAH4 -> x(3)
% EDTAH2 -> x(4)
% EDTAH  -> x(5)
% EDTA   -> x(6)
% Fe     -> x(7)
% FeEDTA -> x(8)
% H      -> x(9)
%----------------

Ksp  = 2.64e-38; 
k1   = 10^-2;
k2   = 10^-2.67;
k3   = 10^-6.16;
k4   = 10^-10.26;
k    = 10^-25;
EDTA = 0.1;
pH   = 13;

y = zeros(size(x));

y(1) = k1*x(3) - x(1)*x(9);
y(2) = k2*x(1) - x(4)*x(9);
y(3) = k3*x(4) - x(5)*x(9);
y(4) = k4*x(5) - x(6)*x(9);
y(5) = k *x(8) - x(7)*x(6);
y(6) = EDTA-x(8)-x(6)-x(5)-x(4)-x(3)-x(1);
y(7) = Ksp - x(7)*x(2)^3;
y(8) = x(9) - 10^-pH;
y(9) = x(2) - 10^-(14-pH); 
