How can interpret the sensitivity analysis of an ordinary differential equation?

Question

Suppose we have the following system of differential equations:

$$ \dot{x_1} = - c_1x_1x_2 + c_2x_3 -c_3x_1 +c_4x_2 $$

$$ \dot{x_2} = - c_1x_1x_2 + c_2x_3 + c_3x_1 -c_4x_2 $$

$$ \dot{x_3} = c_1x_1x_2 - c_2x_3 $$ where all $c_i$ are positive. Suppose we want to compute sensitivity of trajectories with respect to $c_3$. Let's define $\xi_i=\frac{\partial x_i(t)}{\partial c_3}$. Then we have the following equations:

$$ \dot{\xi_1} = -c_1\xi_1x_2 -c_1x_1\xi_2 + c_2\xi_3 -x_1 -c_3 \xi_1 + c_4 \xi_2 $$

$$ \dot{\xi_2} = -c_1\xi_1x_2 -c_1x_1\xi_2 + c_2\xi_3 + x_1 + c_3 \xi_1 - c_4 \xi_2 $$

$$ \dot{\xi_3} = c_1 \xi_1x_2 + c_1 x_1 \xi_2 -c_2 \xi_3 $$

My question is why the perturbation around $c_3 = 3 \pm 0.2$ using uniform random variable is not matched with the solution of the system of equations and sensitivities. To explore that, I have done the following:

If we solve the above 6 differential equations simultaneously using the following codes including the function and values for $c_i$:

function dydt = conc(t,y)
c_1 = 1;
c_2 = 2;
c_3 = 3;
c_4 = 4;

dydt = [-c_1*y(1)*y(2) + c_2*y(3) - c_3*y(1) + c_4*y(2);...
    -c_1*y(1)*y(2) + c_2*y(3) + c_3*y(1) - c_4*y(2);...
    c_1*y(1)*y(2) - c_2*y(3);...
    -c_1*y(4)*y(2) - c_1*y(1)*y(5) + c_2*y(6) - y(1) - c_3*y(4) + c_4*y(5);...
    -c_1*y(4)*y(2) - c_1*y(1)*y(5) + c_2*y(6) + y(1) + c_3*y(4) - c_4*y(5);...
    c_1*y(4)*y(2) + c_1*y(1)*y(4) - c_2*y(6)];

,

clc
clear
close
[t,y] = ode45(@conc,[0 1],[1;1;1;0;0;0]);
figure (1)
plot(t,y(:,1),'-.',t,y(:,2),'-+',t,y(:,3),'-o')
title({'Solution for $x(t)$ with ODE45'},'Interpreter','latex');
xlabel('Time t');
ylabel({'$x(t) = (x_1(t),x_2(t),x_3(t))$'},'Interpreter','latex');
legend({'$x_1(t)$','$x_2(t)$','$x_3(t))$'},'Interpreter','latex')
figure (2)
plot(t,y(:,4),'-.',t,y(:,5),'-+',t,y(:,6),'-o')
title({'Solution for $\frac{\partial x(t)}{\partial c_3}$ with             
ODE45'},'Interpreter','latex');
xlabel('Time t');
ylabel({'$\frac{\partial x(t)}{\partial c_3}=(\xi_1,\xi_2,\xi_3) 
$'},'Interpreter','latex');
legend({'$\frac{\partial x_1(t)}{\partial c_3}$',...
'$\frac{\partial x_2(t)}{\partial c_3}$',...
'$\frac{\partial x_3(t)}{\partial c_3}$'},'Interpreter','latex')

We have the following responses for initial conditions $x_0=(1,1,1)^{\top}$ and $\xi_0=(0,0,0)^{\top}$.

However, if we plot 10 random trajectories which are obtained by perturbing $c_3$ we have the following curves for $x_i$'s:

function dydt = conc_ext(t,y)
c_1 = 1;
c_2 = 2;
c_3 = 2.8 + .4*rand;
c_4 = 4;

dydt = [-c_1*y(1)*y(2) + c_2*y(3) - c_3*y(1) + c_4*y(2);...
        -c_1*y(1)*y(2) + c_2*y(3) + c_3*y(1) - c_4*y(2);...
        c_1*y(1)*y(2) - c_2*y(3)];

,

clc
clear
close
n = 10; 
for i = 1:n
[t,y] = ode45(@conc_ext,[0 1],ones(3,1));
plot(t,y(:,1),'-.',t,y(:,2),'--',t,y(:,3),'.')
hold on
end

% for i = 1:n
% hold on
% end

title({'Solution for $x(t)$ with ODE45'},'Interpreter','latex');
xlabel('Time t');
ylabel({'$x(t) = (x_1(t),x_2(t),x_3(t))$'},'Interpreter','latex');
legend({'$x_1(t)$','$x_2(t)$','$x_3(t))$'},'Interpreter','latex')

The responses are the following:

As the last figure demonstrates $x_3$ is almost constant as $c_3$ perturbed?

Why this happens?