RMS errors of integration

I am simulating the dynamics of an autonomous system described by x1dot=-ax1+x2, x2dot=-bx2+u, udot=-cu. I have generated 100 random variables for x1, x2, u, a, b, and c and then ran 100 iterations of the model. I then plotted x1(t), x2(t), and u(t) with given initial conditions [0 0 0]=[x1 x2, u]. My code is posted below. I am curious how to get the RMS of the errors for each of the integrations for all of the 100 iterations that make up the plots of x1(t), x2(t), and u(t). Any help on this issue would be greatly appreciated. Thanks.

low=0;%constants a,b,c are all positive constants so this is the initial lower bound for what their values can be

high=1000;%for purpose of simplifying simulation upper bound of values for x1, x2, u, a, b, and c was set to 1000. Theoretically, this bound can be almost infinite and those variables can take on larger values

ASa=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1000 for constant a

lowb=1/a;%ensures b > 1/a so vdot is negative definite

ASb=(high-lowb).*rand(100,1) + lowb;%column vector of 100 uniformly distributed decimals between 0 and 1000 for constant b

lowc=1/(2*b);%ensures c > 1/2b so vdot is negative definite

ASc=(high-lowc).*rand(100,1) + lowc;%column vector of 100 uniformly distributed decimals between 0 and 1000 for constant c

ASx1=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1000 for variable x1

ASx2=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1000 for variable x2

ASu=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1000 for variable u

for ii = 1:100 %signifies 100 iterations for each variable x1,x2,u,a,b,c

a =ASa(ii);%a takes on a specific value from the 100 randomly generated values for each iteration

b=ASb(ii);%b takes on a specific value from the 100 randomly generated values for each iteration

c=ASc(ii);%c takes on a specific value from the 100 randomly generated values for each iteration

x1=ASx1(ii);%x1 takes on a specific value from the 100 randomly generated values for each iteration

x2=ASx2(ii);%x2 takes on a specific value from the 100 randomly generated values for each iteration

u=ASu(ii);%u takes on a specific value from the 100 randomly generated values for each iteration

tspan = [0 5];%time span for elapsed simulation

options = odeset('reltol', 1e-5, 'abstol', 1e-8 );%default options for relative and absolute tolerance. The tolerances are used to limit the local discretization error during integration

x0 = [0 0 0];%initial conditions for x1, x2, and u

[t, y] = ode15s(@(t,x) dynamics(t,x,u,a,b,c), tspan, x0, options );%implementation of ode function in matlab to plot our system dynamics described in the function beneath the end of this for loop

figure(1)

plot(t,y(:,1))%plotting of x1(t)

xlabel('Time')

ylabel('x1(t)')

hold on

figure(2)

plot(t,y(:,2))%plotting of x2(t)

xlabel('Time')

ylabel('x2(t)')

hold on

figure(3)

plot(t,y(:,3))%plotting of u(t)

xlabel('Time')

ylabel('u(t)')

hold on

end

hold off

function system = dynamics(t,x,u,a,b,c)

system = zeros(3,1);%initialize system matrix

system(1) = -a*x(1)+x(2);%first row of system matrix defined by equation for dx1/dt

system(2) = -b*x(2)+u;%second row of system matrix defined by equation for dx2/dt

system(3) = -c*u;%third row of system matrix defined by equation for du/dt

end