function out=ExempleReconNonLineaire()
compteur=0; critere=Inf;

syms d1 d2 d3 x1 x2 x3 lambda1 lambda2
%Lagrangien augmenté
% L(d1,d2,d3,x1,x2,x3,lambda1,lambda2)=(13.5-d1)^2/.5^2+(16.1-d2)^2/.5^2+(33.2-d3)^2/.3^2+(190.35-x1*d1)^2/8.13^2+(341.32-x2*d2)^2/19.28^2+(999.32-x3*d3)^2/199.41^2+lambda1*(d1+d2-d3)+lambda2*(x1*d1+x2*d2-x3*d3);
L(d1,d2,d3,x1,x2,x3,lambda1,lambda2)=(13.5-d1)^2/.5^2+(16.1-d2)^2/.5^2+(33.2-d3)^2/.3^2+(14.1*d1-x1*d1)^2/8.13^2+(21.2*d2-x2*d2)^2/19.28^2+(30.1*d3-x3*d3)^2/199.41^2+lambda1*(d1+d2-d3)+lambda2*(x1*d1+x2*d2-x3*d3);

% les 8 équations venant de la dérivation du langragien augmenté:
dL1(d1,d2,d3,x1,x2,x3,lambda1,lambda2)=diff(L(d1,d2,d3,x1,x2,x3,lambda1,lambda2),d1)
dL2(d1,d2,d3,x1,x2,x3,lambda1,lambda2)=diff(L(d1,d2,d3,x1,x2,x3,lambda1,lambda2),d2)
dL3(d1,d2,d3,x1,x2,x3,lambda1,lambda2)=diff(L(d1,d2,d3,x1,x2,x3,lambda1,lambda2),d3)
dL4(d1,d2,d3,x1,x2,x3,lambda1,lambda2)=diff(L(d1,d2,d3,x1,x2,x3,lambda1,lambda2),x1)
dL5(d1,d2,d3,x1,x2,x3,lambda1,lambda2)=diff(L(d1,d2,d3,x1,x2,x3,lambda1,lambda2),x2)
dL6(d1,d2,d3,x1,x2,x3,lambda1,lambda2)=diff(L(d1,d2,d3,x1,x2,x3,lambda1,lambda2),x3)
dL7(d1,d2,d3,x1,x2,x3,lambda1,lambda2)=diff(L(d1,d2,d3,x1,x2,x3,lambda1,lambda2),lambda1)
dL8(d1,d2,d3,x1,x2,x3,lambda1,lambda2)=diff(L(d1,d2,d3,x1,x2,x3,lambda1,lambda2),lambda2)

%matrice jacobienne
DF(d1,d2,d3,x1,x2,x3,lambda1,lambda2)=jacobian([dL1, dL2, dL3, dL4,dL5, dL6, dL7, dL8],[d1,d2,d3,x1,x2,x3,lambda1,lambda2]);
%vecteur contenant les équations
FF(d1,d2,d3,x1,x2,x3,lambda1,lambda2)=[dL1(d1,d2,d3,x1,x2,x3,lambda1,lambda2),dL2(d1,d2,d3,x1,x2,x3,lambda1,lambda2),dL3(d1,d2,d3,x1,x2,x3,lambda1,lambda2),dL4(d1,d2,d3,x1,x2,x3,lambda1,lambda2),dL5(d1,d2,d3,x1,x2,x3,lambda1,lambda2),dL6(d1,d2,d3,x1,x2,x3,lambda1,lambda2),dL7(d1,d2,d3,x1,x2,x3,lambda1,lambda2),dL8(d1,d2,d3,x1,x2,x3,lambda1,lambda2)]';

% estimé initial (mesures des bilans non-ajustés)
d1=13.5; d2=16.1; d3=33.2; x1=14.1; x2=21.2; x3=30.1; lambda1=10;lambda2=.01;

%méthode de Newton matricielle
while critere> 1e-10;
    sol=[d1,d2,d3,x1,x2,x3,lambda1,lambda2]'  %vecteur solution
    solnew=sol-eval(DF(d1,d2,d3,x1,x2,x3,lambda1,lambda2))\eval(FF(d1,d2,d3,x1,x2,x3,lambda1,lambda2)); % nouveau vecteur solution
    
    critere=norm(sol-solnew,2)/norm(sol,2); % epsilon
    d1=solnew(1); d2=solnew(2); d3=solnew(3); x1=solnew(4); x2=solnew(5); x3=solnew(6); lambda1=solnew(7);lambda2=solnew(8);
    compteur=compteur+1
    epsilon(compteur)=critere;
end
sol=[d1,d2,d3,x1,x2,x3,lambda1,lambda2]'
compteur
% impression de la convergence de la méthode de Newton
hf1=figure(1); 
semilogy(epsilon,'bx-','LineWidth',1.5)
ylabel('{\epsilon} [-]','FontSize',20)
xlabel('Itérations [-]','FontSize',20)
set(gca,'LineWidth',2,'FontSize',18)
title('Convergence de la méthode de Newton')
box on;
axis([1 compteur 1e-15 1e4 ])
text(1.2,10^(-7),['D_{1,a}=',num2str(sprintf('%.3f',sol(1))),' t/h'],'FontSize',18)
text(1.2,10^(-9),['D_{2,a}=',num2str(sprintf('%.3f',sol(2))),' t/h'],'FontSize',18)
text(1.2,10^(-11),['D_{3,a}=',num2str(sprintf('%.3f',sol(3))),' t/h'],'FontSize',18)
text(1.2,10^(-13),['{\lambda_1}=',num2str(sprintf('%.3f',sol(7))),' h/t'],'FontSize',18)
text(3.2,10^(2.5),['x_{1,a}=',num2str(sprintf('%.3f',sol(4))),' %'],'FontSize',18)
text(3.2,10^(0.5),['x_{2,a}=',num2str(sprintf('%.3f',sol(5))),' %'],'FontSize',18)
text(3.2,10^(-1.5),['x_{3,a}=',num2str(sprintf('%.3f',sol(6))),' %'],'FontSize',18)
text(3.2,10^(-3.5),['{\lambda_2}=',num2str(sprintf('%.3f',sol(8))),' h/t/%'],'FontSize',18)
end