function u = transport(c, T, xmin,xmax, dx, dt)
%%%%%%%%%%%%%%%%%%%
%ce script résout l'équation de transport linéaire
% u_t + c u_x = 0
% u(x,O) = f(x)
%
% sur  xmin<x<xmax, 0<t<T
% avec la méthode des différences finies
% on code 4  schémas explicites:
% discrétisation amont
% discrétisation centrée
% Lax-Friedrichs
% Lax-Wendroff 
% on affiche la CFL= c*dt/dx
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% maillage
x=(xmin:dx:xmax);
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);
u(:,1,4) = f(x);
% CFL
r= c*dt/dx
% boucle en temps
for n=1:length(t)
    %%%%%%%%schema  décentré amont%%%%%%%%
%     %affectation vectorielle plus rapide qu'une boucle
     j=(2:length(x));
     u(j,n+1,1) = r*u(j-1,n,1) + (1-r)* u(j,n,1);
%     %condition au bord x=xmin periodique
%    %u(1,n+1,1) = r*u(end,n,1) + (1-r)* u(1,n,1);
%    % condition au bord x = xmin constante
     u(1,n+1,1) = u(1,n,1);
%     %%%%%%%%schema centré%%%%%%%%%
     j=(2:length(x)-1);
     u(j,n+1,2) = (r/2)*u(j-1,n,2)+ u(j,n,2)...
                                -(r/2)*u(j+1,n,2);
%     %condition au bord x=xmin periodique
%     u(1,n+1,2) = (r/2)*u(end,n,2)+ u(1,n,2) - (r/2)*u(2,n,2);
%     %condition au bord x=xmax periodique
%     u(end,n+1,2) = (r/2)*u(end-1,n,2)+ u(end,n) - (r/2)*u(1,n,2);
%     % conditions au bord constantes
        u(1,n+1,2) = u(1,n,2);
        u(end,n+1,2) = u(end,n,2);
        %%%%%%%schéma Lax-Friedrichs%%%%%%%
        j=(2:length(x)-1);
       u(j,n+1,3) = (1+r)/2*u(j-1,n,3)+ (1-r)/2* u(j+1,n,3);
%         %condition au bord x=xmin periodique
%         u(1,n+1,3) = (1+r)/2*u(end,n,3)+ (1-r)/2*u(2,n,3);
%         %condition au bord x=xmax periodique
%         u(end,n+1,3) = (1+r)/2*u(end-1,n,3)+(1-r)/2*u(1,n,3);
%        % conditions au bord constantes
          u(1,n+1,3) = u(1,n,3);
          u(end,n+1,3) = u(end,n,3);
%       %%%%%%%schéma Lax-Wendroff%%%%%%%%
        j=(2:length(x)-1);
        u(j,n+1,4) = (r^2+r)/2*u(j-1,n,4)+...
                (1-r^2)*u(j,n,4) + (r^2-r)/2* u(j+1,n,4);
%             %condition au bord x=xmin periodique
%             u(1,n+1,4) = (r^2+r)/2*u(end,n,4)+ ...
%                 (1-r^2)*u(1,n,4) + (r^2-r)/2* u(2,n,4);
%             %condition au bord x=xmax periodique
%             u(end,n+1,4) = (r^2+r)/2*u(end-1,n,4)+ ...
%                 (1-r^2)*u(end,n,4) + (r^2-r)/2* u(1,n,4);     
   %conditions au bord constantes
           u(1,n+1,4) = u(1,n,4);
           u(end,n+1,4) = u(end,n,4);
    
end
% tracé des courbes pour chaque schéma
close all;
figure;
subplot(221)
xlim([xmin xmax])
ylim([-0.5 1.5])
hold on
plot(x,f(x),'-b',x, u(:,end,1),'-k');
plot(x, f(x-c*T),'--k');
legend('t=0','t=T','sol. exacte t=T','location','SW');
xlabel('x');
ylabel('u(.,t)');
title('schema upwind',['CFL=' num2str(r)])
%
subplot(222)
xlim([xmin xmax])
ylim([-0.5 1.5])
hold on
plot(x,f(x),'-b',x, u(:,end,2),'-k');
plot(x, f(x-c*T),'--k');
legend('t=0','t=T','sol. exacte t=T','location','SW');
xlabel('x');
ylabel('u(.,t)');
title('schema centré ',['CFL=' num2str(r)])
%
subplot(223)
xlim([xmin xmax])
ylim([-0.5 1.5])
hold on
plot(x,f(x),'-b',x, u(:,end,3),'-k');
plot(x, f(x-c*T),'--k');
legend('t=0','t=T','sol. exacte t=T','location','SW');
xlabel('x');
ylabel('u(.,t)');
title('schéma Lax-Friedrichs ',['CFL=' num2str(r)])
%
subplot(224)
xlim([xmin xmax])
ylim([-0.5 1.5])
hold on
plot(x,f(x),'-b',x, u(:,end,4),'-k');
plot(x, f(x-c*T),'--k');
legend('t=0','t=T','sol. exacte t=T','location','SW');
xlabel('x');
ylabel('u(.,t)');
title('schéma Lax-Wendroff ',['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_(0,1/2)
    %y = 0.*(x<0) + 1.*(0<=x).*(x<1/2);
    % cloche de Gauss.
    % y = exp(-25*( x- (xmin + xmax)/2).^2);
    %
    
end
end