//
// equilibre d'un système chimique (3 especes, 2 composants, 4000 degrés)
//
function y=phi(x)
  g1=-21.89;
  g2=-31.44;
  g3=-38.17;
  y=x;
  X1=x(1);X2=x(2);X3=x(3);l1=x(5);l2=x(6);N=x(4);
  y(1)=log(X1)+g1-2*l1;
  y(2)=log(X2)+g2-2*l2;
  y(3)=log(X3)+g3-2*l1-l2;
  y(4)=X1+X2+X3-1;
  y(5)=2*X1+2*X3-4/N
  y(6)=2*X2+X3-2/N    
endfunction
function M=Jac(x)
  M=zeros(6,6);
  X1=x(1);X2=x(2);X3=x(3);l1=x(5);l2=x(6);N=x(4);
  M(1,1)=1/X1;M(1,5)=-2;
  M(2,2)=1/X2;
  M(2,6)=-2;
  M(3,3)=1/X3;
  M(3,5)=-2;
  M(3,6)=-1;
  M(4,1)=1;
  M(4,2)=1;
  M(4,3)=1;
  M(5,1)=2;
  M(5,3)=2;
  M(5,4)=4/N^2;
  M(6,2)=2;
  M(6,3)=1;
  M(6,4)=2/N^2;
endfunction
Z=[0.3;0.3;0.4;3;0;0];
for i=1:10
    Z=Z-inv(Jac(Z))*phi(Z)
    disp(cond(Jac(Z)));disp(Z)
end
disp('état d''equilibre');disp(Z(1:4))
