function U = calcRefroidissement(T,u0,r,Uext,methode)
% function U = calcRefroidissement(T,u0,r,Uext,methode)
%-----------------------------------------------------------------------------
% calcule la temperature u(t) [°C] du café dans une tasse, obéissant à
% l'équation différentielle ordinaire:     du/dt = -r*[u(t)-uext(t)]   pour t>=0
% avec la condition initiale:              u(0)=u0.
% Le coefficient de la conductivité thermique de la paroi de tasse, r [1/s],
% donne le taux d'échange de l'énergie thermique entre la tasse du café et l'extérieur.
% La température extérieure uext(t) est donnée (et pas nécessairement constante).
% On calcule u(t) résultant aux instants T = [t0,t1,t2,t3,...tn]', où t0=0
% et t1,t2,... n'est pas nécessairement équidistant
% La température extérieure uext(t) est donnée par ses valeurs aux
% instants t ∈ T, données par le vecteur Uext = [uext(t0),uext(t1),...,uext(tn)]';
% La résolution numérique de l'équation différentielle utilise une des
% méthodes (methode) suivantes:
%   "EulerExplicite": méthode d'Euler explicite (progressive) d'ordre
%      de précision 1 en pas de temps dt,
%   "Heun": méthode de Heun (d'Euler modifiée) d'ordre de orécision 2,
%   "RungeKutta4": est la méthode de Runge-Kutta 4 (classique)
%   "EulerImplicte": méthode d'Euler implicite (rétrograde) d'ordre 1 (la
%      seule méthode implicite parmi les méthodes présentes)
%*** paramètres:
% T(n,:)        -in-    instants de temps T=[t0,t1,t2,...,tn]'
% u0            -in-    température initiale du café en [°C]
% r             -in-    coefficient d'échange thermique [1/s]
% Uext(n,:)     -in-    Uext=[uext(t0),uext(t1),...]
% methode       -in-    string avec la méthode de résolution de l'équation diff.
% U(n,:)        -out-   la température du café aux instants T=[t0,t1,...,tn]'
%-----------------------------------------------------------------------------

%*** init:
  U = 0*T;
  U(1) = u0;
  n = length(T);
%*** uext(t) est obtenue en interpolant les valeurs (T,Uext), si besoin:
  uext = @(t) interp1(T,Uext,t,'pchip');
%*** la forme générale d'un problème de Cauchy est du/dt = f(t,u) avec u(t0)=u0:
  f = @(t,u) -r*(u-uext(t));
%=== EXPRIMEZ f(t,u) ici ^^^^^^

%=== PROGRAMMEZ ICI le méthodes de résolution de l'équation différentielle:
  switch (methode)
  case "EulerExplicite",
    for k=1:n-1
      dt = T(k+1)-T(k);
      U(k+1) = U(k)+dt*f(T(k),U(k));
    end;
  case "Heun",
    for k=1:n-1
      dt = T(k+1)-T(k);
      p1 = f(T(k),U(k));
      p2 = f(T(k+1),U(k)+dt*p1);
      U(k+1) = U(k) + (dt/2)*(p1+p2);
    end;
    case "RungeKutta4",
    for k=1:n-1
      dt = T(k+1)-T(k);
      p1 = f(T(k),U(k));
      p2 = f(T(k) + dt/2, U(k)+(dt/2)*p1);
      p3 = f(T(k) + dt/2, U(k)+(dt/2)*p2);
      p4 = f(T(k) + dt, U(k)+ dt*p3);
      U(k+1) = U(k) + (dt/6)*(p1+2*p2+2*p3+p4);
    end;
  case "EulerImplicite",
    for k=1:n-1
      dt = T(k+1)-T(k);
      U(k+1) = (U(k)+dt*r*uext(T(k+1)))/(1+dt*r);
    end;
  otherwise,
    error('cette méthode n´est pas implémentée');
  end;
%=== FIN DE SOLUTION
  return;
end