function Y=resout_euler_implicite(N,T,y0,fcn,nmax,epsilon)
% 
%
%
%	resoud_euler_implicite : résolution d'une équation différentielle par la méthode d'Euler implicite.
%
% *********************************************************
%
%    X=resout_euler_implicite(N,T,y0,fcn,nmax,epsilon) calcule les valeurs aux instants 0,h,2h,...,Nh=T par la méthode
%    d'Euler implicite pour l'équation différentielle y't)=f(t,y(t)) sur [0,T] et y(0)=y0.
%    A chaque itération, l'équation y_{n+1}=y_n+hf(t_{n+1},y_{n+1}) est résolue par la méthode de Newton
%    où la valeur initiale est choisie égale à y_n.
%    
%    La dérivée partielle df/dy est calculée en symbolique.
%
%	* variables d'entrées :
%      N : nombre de valeurs calculées ; 
%      T : borne supérieure de l'intervalle ;
%      y0 : valeur de y à t=0 ;
%      fcn : chaine de caractères (de type 'f(t,y)') représentant la fonction f
%      epsilon : précision souhaitée (pour la méthode de Newton)
%      nmax : nombre maximal d'itérations (pour la méthode de Newton)
%
%	* variables de sortie :
%     Y : valeurs calculée aux instants 0,h,2h,...,Nh=T.
%
% 
% 
% ************ Fonctions auxiliaires utilisées ************
%
%       aucune
%
% *********************************************************
%

% Contrôles d'entrée

% nombre d'arguments
if nargin~=6
   error('nombre d''arguments de la fonction incorrect');
end
   
% Corps d'algorithme
h=T/N;
Y=zeros(1,N+1);
Y(1)=y0;
auxiY=y0;
syms t y;
g=eval(fcn);
fcnpsymb=diff(g,y);
fcnp=char(fcnpsymb);
yapp=y0;
for i=1:N
   t=h*i;
   y=yapp;
   fx=-h*eval(fcn);
   fxp=1-h*eval(fcnp);
   test=((fx)~=0);
   n=0;
   v=yapp;
   while(test)
      oldv=v;
      v=v-fx/fxp;
      y=v;
      fx=v-yapp-h*eval(fcn);
      fxp=1-h*eval(fcnp);
      n=n+1;
      test=~((fx==0)|(n>=nmax)|(abs(v-oldv)<=epsilon));
   end   
   yapp=v;
   if (n>=nmax)
      error('attention, n>nmax pour la méthode de Newton : arrêt du programme');
   end
   Y(i+1)=yapp;
end




