function u = diffusion(c,T, xmin,xmax, dx, dt)
%%%%%%%%%%%%%%%%%%%
%ce script résout l'équation de diffusion
% u_t - c u_xx = 0
% u(x,O) = f(x)
%
% sur  xmin<x<xmax, 0<t<T
% avec la méthode des différences finies
% on code 3  schémas
% discrétisation explicite centrée
% discrétisation implicite centrée
% discrétisation de Crank-Nicolson
% 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 3 schémas s=1 explicite,
%2 implicite ,3 Crank-Nicolson
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% maillage
x=(xmin:dx:xmax);
N=length(x);
t=(0:dt:T);
%initialisation pour les 4 schémas
u=zeros(length(x),length(t),4);
% donnee initiale pour les 4 schémas.
u(:,1,1) = f(x);
u(:,1,2) = f(x);
u(:,1,3) = f(x);
% CFL
r= c*dt/dx^2
%%%%%%%%schema  explicite %%%%%%%%
% boucle en temps
for n=1:length(t)
    %affectation vectorielle plus rapide
    j=(2:N-1);
    u(j,n+1,1) = r*u(j-1,n,1) + ...
        (1-2*r)* 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)* 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)* 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)*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);
    u(:,n+1,2) = A\u(:,n,2);
end
%%%%%%%schéma Crank-Nicolson%%%%%%%
%%%création des  matrices tridiagonales
%%%%% au format sparse
A = spdiags([-r/2*e (1+r)*e -r/2*e],-1:1,N,N);
%modification des coeff à cause des conditions
%limites périodiques
A(1,N)= -r/2;
A(N,1)= -r/2;
%%%%%%%%%%%%%%%%%%%%%
B = spdiags([r/2*e (1-r)*e r/2*e],-1:1,N,N);
%modification des coeff à cause des conditions
%limites périodiques
B(1,N)= r/2;
B(N,1)= r/2;
%boucle en temps
% a chaque itération on résout A U^{n+1} = B U^n
for n=1:length(t);
    u(:,n+1,3) = A\(B*u(:,n,3));
end
% tracé des courbes pour chaque schéma
close all;
figure;
subplot(131)
xlim([xmin xmax])
ylim([-0.5 1.5])
hold on
plot(x,f(x),'-b',x, u(:,end,1),'-k');
legend('t=0','t=T','location','SW');
xlabel('x');
ylabel('u(.,t)');
title('schéma explicite',['CFL=' num2str(r)])
% %
subplot(132)
xlim([xmin xmax])
ylim([-0.5 1.5])
hold on
plot(x,f(x),'-b',x, u(:,end,2),'-k');
legend('t=0','t=T','location','SW');
xlabel('x');
ylabel('u(.,t)');
title('schéma implicite ',['CFL=' num2str(r)])
% %
subplot(133)
xlim([xmin xmax])
ylim([-0.5 1.5])
hold on
plot(x,f(x),'-b',x, u(:,end,3),'-k');
legend('t=0','t=T','location','SW');
xlabel('x');
ylabel('u(.,t)');
title('schéma Crank-Nicolson ',['CFL=' num2str(r)])
%
% 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);
        %
        
    end
end