function u = difrea(a, c, T, xmin,xmax, dx, dt)
%%%%%%%%%%%%%%%%%%%
%ce script résout l'équation de diffusion réaction
% u_t - c u_xx - a u = 0
% u(x,O) = f(x)
%
% sur  xmin<x<xmax, 0<t<T
% avec la méthode des différences finies
% on code 2  schémas
% discrétisation explicite centrée
% discrétisation implicite centrée
%
% on affiche la CFL parabolique = c*dt/dx^2
%
%le script retourne la matrice u(j,n,s) = u(x_j, t_n,s)
% pour les 2 schémas s=1 explicite,
% s=2 implicite 
% dans le cas où f(x)=sin(x) sur (-pi,pi) on peut comparer à 
% la solution exacte sin(x)*exp(a-c)*t
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% maillage
x=(xmin:dx:xmax);
N=length(x);
t=(0:dt:T);
%initialisation pour les 2 schémas
u=zeros(length(x),length(t),2);
% donnee initiale pour les 2 schémas.
u(:,1,1) = f(x);
u(:,1,2) = f(x);
% CFL
r= c*dt/dx^2;
s=a*dt;
%%%%%%%%schema  explicite %%%%%%%%
% boucle en temps
for n=1:length(t)-1
    %affectation vectorielle plus rapide
    j=(2:N-1);
    u(j,n+1,1) = r*u(j-1,n,1) + ...
        (1-2*r+s)* u(j,n,1) + r*u(j+1,n,1);
    %condition au bord x=xmin periodique
    u(1,n+1,1) = r*u(end,n,1) +  ...
        (1-2*r +s)* u(1,n,1) + r*u(2,n,1);
    %condition au bord x=xmax periodique
    u(end,n+1,1) = r*u(end-1,n,1) + ...
        (1-2*r + s)* u(end,n,1) + r*u(1,n,1);
end
%%%%%%%%%schema implicite%%%%%%%%%
e = ones(N,1);
%%%création de la matrice tridiagonale
% on utilise le format sparse (creux)
A = spdiags([-r*e (1+2*r-s)*e -r*e],-1:1,N,N);
%modification des coeff à cause des conditions
%limites périodiques
A(1,N)= -r;
A(N,1)= -r;
% boucle en temps
% a chaque itération on résout A U^{n+1} = U^n
for n=1:length(t)-1;
    u(:,n+1,2) = A\(u(:,n,2));
end

% tracé des courbes pour chaque schéma
%close all;
figure();
subplot(121)
xlim([xmin xmax])
%ylim([-2 2])
hold on
plot(x,f(x),'-b',x, u(:,end,1),'-k',x, exp((a-c)*T).*sin(x),'or');
%legend('t=0','t=T','location','SW');
legend('t=0','t=T','exacte');
xlabel('x');
ylabel('u(.,t)');
title("schema explicite")
%title('schéma explicite',['CFL=' num2str(r)])
% %
subplot(122)
xlim([xmin xmax])
%ylim([-2 2])
hold on
plot(x,f(x),'-b',x, u(:,end,2),'-k',x, exp((a-c)*T).*sin(x),'or');
%legend('t=0','t=T','location','SW');
legend('t=0','t=T','exacte');
xlabel('x');
ylabel('u(.,t)');
%title('schéma implicite ',['CFL=' num2str(r)])
title("schema implicite")
%
%
% donnée initiale
    function y=f(x)
        %decommenter la fonction choisie.
        % fonction Heaviside lissée
        % y = 0.*(x<0)+4*x.^2.*(3-4*x).*(0<=x).*(x<=0.5)+ 1.*(x>0.5);
        %fonction 1_(x<0)
        % y = 1.*(x<0);
        %fonction creneau 1_(-1/2,1/2)
        % y =  1.*(-0.5<=x).*(x<1/2);
        % cloche de Gauss.
        % y = exp(-25*( x- (xmin + xmax)/2).^2);
        %
        % sinus
          y= sin(x);
        % fourier
        %  y= sin(x) + 0.5*sin(2*x);
        
    end
end