// Schéma implicite 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) = 1
// c(t,L) = 0


// Paramètres
w = 0.9;
v = 0.003;
u = 0.1;
s0 = 2;
k = 10000;
p = 3;
CBg = 1;
sig = u/(w+s0/p);

// Discrétisation 
M = 100;
L = 1;
dx = L/M;
x = (0:M)'*dx;
T = 15;
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;

// Méthode de Newton
function Z = Newton(C,S,Lap,Diff,E,dt,k,p,w)
    Z = C;
    M = length(Z);
    for n=1:5
        F0 = -Z+C+dt*(Lap*C-Diff*C-k*(Z.*S)./(1+k*p*dt*Z)+E)/w;
        DF0 = -sparse(eye(M,M))-(k*dt/w)*sparse(diag(S./(1+k*p*dt*Z)-k*p*dt*(Z.*S)./((1+k*p*dt*Z).^2)));
        Z = Z-DF0\F0;
    end
endfunction

// résolution de l'EDP
for n=1:N
    c = Newton(c,s,v*Lap,u*Diff,E,dt,k,p,w);
    s = s./(1+dt*p*k*c);
    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

















