// Schéma explicite pour résoudre
// w c_t = v c_xx - u c_x -kcs
//   s_t = -pkcs
// c(0,x) = 0
// s(0,x) = s0
// c(t,0) = cd(t)
// c(t,L) = 0


// Paramètres
w = 0.9;
v = 0.01;
u = 0.4;
s0 = 2;
k = 1;
p = 10;


// Discrétisation 
M = 200;
L = 4;
dx = L/M;
x = (0:M)'*dx;
T = 10;
dt = 0.01;
N = T/dt;

// Matrice du Laplacien 
Lap = sparse(-2*eye(M-1,M-1)+diag(ones(M-2,1),1)+diag(ones(M-2,1),-1));
Lap = Lap/(dx^2);

// Matrice de transport
Diff = sparse(eye(M-1,M-1)-diag(ones(M-2,1),-1));
Diff = Diff/dx;

// Vecteur condition de bord à gauche
E = [v/(dx^2)+u/dx;zeros(M-2,1)];

// Conditions initiales
CIc = zeros(M-1,1);
CIs = s0*ones(M-1,1);

c = CIc;
s = CIs;
sg = s0;
t = 0;

// résolution de l'EDP
for n=1:N
    CBg = t*exp(-t/2);
    c = (c+dt*(v*Lap*c-u*Diff*c-k*c.*s+CBg*E)/w);
    s = s-dt*p*k*c.*s;
    sg = sg*(1-dt*p*k*CBg);
    drawlater();
    clf();
    a=gca(); 
    a.data_bounds=[0,0;L,s0];
    plot(x,[CBg;c;0],'b',x,[sg;s;s(M-1)],'k','thickness',3);
    xlabel("x")
    drawnow();
    t = t+dt;
end

















