Modèle Gillespie - Plot - MatLab

[Eindhov]

Bonjour tout le monde,
J'ai un projet à rendre sur le modèle Gillespie. Si vous ne connaissez pas de nombreux articles sont présents sur le net.
Mon prblème est que je dois développer ça avec MatLab, outil que j'ai utilisé deux fois dans ma vie...
Dans mon programme j'essaie de plot a en fonction du temps, mais cela ne me retourne rien. Quelqu'un peut)-il m'aider svp ?
Voici le code :

clear all;

% Set initial state
a = 1; 
b = 1;
c = 0;
d = 0;
e = 0;
t = 0;
tend = 1000;

%Mesoscopic : 
r1 = 4;     %A->ARNA
r2 = 1;     %ARNA->ARNA + A
r3 = 1;     %ARNA->
r4 = 0.1;   %A->
r5 = 100;   %A + R->AR (AR->A + R)
r6 = 0.01;  %R->
r7 = 0.1;   %RRNA->R
r8 = 0.02;  %RRNA->
r9 = 0.001; %A->RRNA

while t<tend 
    u1 = r1*a;
    u2 = r2*c;
    u3 = r3*c;
    u4 = r4*a;
    u5 = r5*a*b;
    u6 = r6*b;
    u7 = r7*d;
    u8 = r8*d;
    u9 = r9*a;
    
    lambda = u1+u2+u3+u4+u5+u6+u7+u8+u9;
    Y = rand;
    tau = -log(Y)./lambda;
    t=t+tau;
    
    x1 = u1/lambda;
    x2 = (u1+u2)/lambda;
    x3 = (u1+u2+u3)/lambda;
    x4 = (u1+u2+u3+u4)/lambda;
    x5 = (u1+u2+u3+u4+u5)/lambda;
    x6 = (u1+u2+u3+u4+u5+u6)/lambda;
    x7 = (u1+u2+u3+u4+u5+u6+u7)/lambda;
    x8 = (u1+u2+u3+u4+u5+u6+u7+u8)/lambda;
    
    z = rand;
    
    if z<x1
        a = a - 1;
        c = c +1;
        
    elseif z<x2
        a = a + 1;
       
    elseif z<x3
        c = c - 1;
        
    elseif z<x4
        a = a - 1;
       
    elseif z<x5
        a = a - 1;
        b = b - 1;
        c = c + 1;
        
    elseif z<x6
        b = b-1;
       
    elseif z<x7
        b = b + 1;
        d = d - 1;
        
    elseif z<x8
        d = d - 1;
        
    else 
        a = a - 1;
        d = d-1;
      
    end
end

plot(t,a);
xlabel('Time');
ylabel('Number of particles A');
plot(b,t);

Merci d'avance .
PS : ne tenez pas compte des commentaires à côté des déclarations des variables r .
-----