Issues plotting ode in a for loop

I am simulating the dynamics of an autonomous system which is made up of following equations: x1dot=-ax1+x2, x2dot=-bx2+u, udot=-cu. The variables a,b,c are constants with applied constraints on them, the initial conditions for [x1 x2 u]=[0 0 0]. My code is posted below. I am trying to plot x1(t), x2(t), and u(t) on three separate plots for a given time span and for 100 random iterations of x1,x2,u,a,b,c. The issue I am running into is that my plots keep generating over and over for each iteration instead of placing the 100 trajectories for x1(t), x2(t), and u(t) on three separate plots in which all of the iterations for the given time span are on a single plot. Any help on this would be much appreciated. Thanks.

low=1;
high=1000;
ASa=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for constant a 
ASb=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for constant b
ASc=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for constant c
ASx1=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for variable x1
ASx2=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for variable x2
ASu=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for variable u
for ii = 1:100 %signifies 100 iterations for each variable x1,x2,u,a,b,c
    a =ASa(ii);
    b=ASb(ii); 
    b = max([b, 1/a]); % Make sure b > 1/a
    c=ASc(ii);
    c = max([c, 1/(2*b)]); % Make sure c is > 1/(2b)
    x1=ASx1(ii);
    x2=ASx2(ii);
    u=ASu(ii);
    tspan = [0 5];
    options = odeset('reltol', 1e-5, 'abstol', 1e-8 );
    x0 = [0 0 0];
    [t, y] = ode15s(@(t,x) dynamics(t,x,u,a,b,c), tspan, x0, options );
    figure(1)
    plot(t,y(:,1))
    hold on
    figure(2)
    plot(t,y(:,2))
    figure(3)
    plot(t,y(:,3))
end
title('Problem 1 System Dynamics')
xlabel('Time')
ylabel('Dynamic System')
hold off

function system = dynamics(t,x,u,a,b,c)
system = zeros(3,1);
system(1) = -a*x(1)+x(2);
system(2) = -b*x(2)+u;
system(3) = -c*u;
% system=[-a*x1+x2;-b*x2+u;-c*u] 
end

Answer 1

[Sulaymon Eshkabilov]

There are a couple of "hold on" s are missing and thus, not all simulated data get plotted. Here is the corrected part of the code:

low=1;
high=1000;
ASa=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for constant a 
ASb=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for constant b
ASc=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for constant c
ASx1=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for variable x1
ASx2=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for variable x2
ASu=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for variable u
for ii = 1:100 %signifies 100 iterations for each variable x1,x2,u,a,b,c
    a =ASa(ii);
    b=ASb(ii); 
    b = max([b, 1/a]); % Make sure b > 1/a
    c=ASc(ii);
    c = max([c, 1/(2*b)]); % Make sure c is > 1/(2b)
    x1=ASx1(ii);
    x2=ASx2(ii);
    u=ASu(ii);
    tspan = [0 5];
    options = odeset('reltol', 1e-5, 'abstol', 1e-8 );
    x0 = [0 0 0];
    [t, y] = ode15s(@(t,x) dynamics(t,x,u,a,b,c), tspan, x0, options );
    figure(1)
    plot(t,y(:,1))
    hold on
    figure(2)
    plot(t,y(:,2))
    hold on         % Necessary
    figure(3)
    plot(t,y(:,3))
    hold on         % Necessary
end
title('Problem 1 System Dynamics')
xlabel('Time')
ylabel('Dynamic System')
hold off

function system = dynamics(t,x,u,a,b,c)
system = zeros(3,1);
system(1) = -a*x(1)+x(2);
system(2) = -b*x(2)+u;
system(3) = -c*u;
% system=[-a*x1+x2;-b*x2+u;-c*u] 
end

You may also consider to compute 1st and then plot the data, e.g.:

...
for ii = 1:100 %signifies 100 iterations for each variable x1,x2,u,a,b,c
   ....
    [t, y] = ode15s(@(t,x) dynamics(t,x,u,a,b,c), tspan, x0, options );
    T{ii}=t;
    Y{ii}=y;
end
...
    

Then plot the results:

for ii=1:100
    figure(1)
    plot(T{ii},Y{ii}(:,1))
    hold on
    figure(2)
    plot(T{ii},Y{ii}(:,2))
    hold on
    figure(3)
    plot(T{ii},Y{ii}(:,3))
    hold on
end
...