function [u] = eqdiff2(N = 10, M = 10, h = 5E-3)

%
%  Schema d'Euler retrograde pour un systeme differentiel
%
%  parametres
%
%  N        : nombre d'inconnues du systeme differentiel
%  M        : nombre de pas de temps
%  h        : pas de temps
%  u        : N-vecteur, approximation de la solution du systeme differentiel
%  diag     : coefficient diagonal de la matrice I+hA
%  sdiag    : coefficient sous-diagonal de la matrice I+hA
%  a        : N-vecteur, diagonale de L tq I+hA=LL^T
%  c        : (N-1)-vecteur, sous-diagonale de L tq I+hA=LL^T
%

alpha=(N+1)^2;

%
% condition initiale
%
for i=1:N
  u(i)=1.;
end

%
% decomposition de cholesky de la matrice I+hA : I+hA = LL^T
%
diag = ?;
sdiag = ?;
a(1) = sqrt(?);
for i = 1 : N - 1
  c(i) = ? / a(i);
  a(i+1) = sqrt(? - c(i) * c(i));
end
%
% schema d'Euler retrograde : (I+hA) u^{n+1} = u^{n}
%
for n = 1 : M
%
% resolution du systeme lineaire Ly = u^{n} */
%
  u(1) = u(1) / a(1);
  for i = 1 : N - 1
    u(i+1) = (u(i + 1) - ? * u(i)) / ?;
  end
%
% resolution du systeme lineaire L^T u^{n+1} = y */
%
  u(N) = u(N) / a(N);
  for i = N - 1 : -1 : 1
    u(i) = (u(i) - ? * u(i+1)) / ?;
  end
%
% impression de la norme euclidienne de u^{n+1}
%
  norm2 = 0;
  for i = 1 : N
    norm2 = norm2 + u(i) * u(i);
  end
  norm2
end

end
