function Y=resout_theta_methode(theta,N,T,y0,fcn,nmax,epsilon)
% 
%
%
%	resout_theta_methode : résolution d'une équation différentielle par la theta-méthode.
%
% *********************************************************
%
%    Y=resout_theta_methode(theta,N,T,y0,fcn,nmax,epsilon) calcule les valeurs aux instants 0,h,2h,...,Nh=T par 
%    la theta-méthode
%    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+theta hf(t_{n+1},y_{n+1})+(1-theta)hf(t_{n},y_{n})
%    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 :
%      theta est un réel dans [0,1] ;
%      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~=7
   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
   y=yapp;
   t=h*(i-1);
   auxidud=h*(1-theta)*eval(fcn);
   t=h*i;
   fx=-h*theta*eval(fcn)-auxidud;
   fxp=1-h*theta*eval(fcnp);
   test=((fx)~=0);
   n=0;
   v=yapp;
   while(test)
      oldv=v;
      v=v-fx/fxp;
      y=v;
      t=h*i;
      fx=v-yapp-h*theta*eval(fcn)-auxidud;
      fxp=1-h*theta*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
