global pH EnT ;
int Row ;
Row=57;
double HEn[Row];double H2En[Row];double En[Row];
%double HEn1[Row];double H2En1[Row];double En1[Row];

EnT=1;
  for i=1:Row
      pH=i/4-0.25;
    X0=[0;1;0;pH;EnT];
    %  X0=[HEn(i);H2En(i); En(i);pH;EnT];
   options=optimset('Display','off');
    [X,fval]=fsolve('myfunCuO123',X0,options);
     HEn(i)=X(1);H2En(i)=X(2); En(i)=X(3);  %*100/EnT
   %HEn1(i)=X(1);H2En1(i)=X(2); En1(i)=X(3);
  


end