function F=myfunCuO(X)

global ClT EnT;
A=[  16.9068 -4.606  0 0;  %Cu2+
    19.1637 -4.606  1 0;  18.4959 -4.606  2 0;  18.1735  -4.606  3 0; %CuCl(1-3)
    41.0883 -4.606  0 1;  62.0456 -4.606  0 2;  59.9729  -4.606   0   3;   31.4848   -2.303   0   2;    %Cu(En)1-3;  CuOHEn2   
    -73.9179  4.606  0  0;   -0.8263  -2.303  0  0;   -18.0988  0  0  0;   -46.4257  2.303  0  0;   -76.1344  4.606  0  0;  9.1254  -4.606  0  0;
    %CuO2(2-)                    CuOH+                      Cu(OH)2(aq)       Cu(OH)3(-)                Cu(OH)4(2-)             Cu2(OH)2(2+)     
    22.8642  -2.303  0  1;   38.6351  -4.606  0  1;  %HEn+     H2En2+
    0 -2.303 0 0;           -32.2348 2.303 0 0]; % H+   OH-

B=[1; X(1); log(X(2)); log(X(3))];  %X(1)=pH      X(2)=Cl      X(3)=En

C=exp(A*B);


D=[1;1;1;1;1;  1;1;1;1;1;  1;1;1;1;0;  0;0;0];

DEn=[0;0;0;0;0;0;0;0;1;0;0];

DH=[0;0;0;0;0;0;0;0;0;1;0];

DOH=[0;0;0;0;0;0;0;0;0;0;1];

DN=[2;1;0;-1;2;  2;2;1;-2;1;  0;-1;-2;2;1; 2;1;-1];


F(1)=dot(D,C)-X(4);

F(2)=X(2)+dot(A(:,3),C)-ClT;

F(3)=X(3)+dot(A(:,4),C)-EnT;

F(4)=dot(DN,C)-X(2);

%X0=[6.0;0.100;0.0500;1.00;0.1;0.1];
%options=optimset('Display','iter');
%[X,fval]=fsolve('mynewfun',X0,options);