Edit: After understanding what you really want in your latest clarification. There is a simpler and intuitive way to code the program such that the simulation only run for approximately 3 or 4 times the Settling Time, 
.
. This allows 
or 
of the plot window to show the transient trajectories of the states from the initial values to the steady-state values. This method is only meaningful for systems that have stable equilibrium points.
or of the plot window to show the transient trajectories of the states from the initial values to the steady-state values. This method is only meaningful for systems that have stable equilibrium points.Analysis shows that your system has a stable equilibrium point at 
and 
and three other unstable equilibrium points at 
, 
, and 
.
and and three other unstable equilibrium points at , , and . Note: I changed the initial values for 
and 
because 
and 
, and to show you the difference between the settling time approach and the event function approach.
and because and , and to show you the difference between the settling time approach and the event function approach.alpha = 0.5;
beta = 0.5;
r1 = 2;
r2 = 3;
s1 = 1;
s2 = 1;
t0 = 0;
tfinal = 5; % Adjust this parameter roughly 4 times the Settling Time, Ts
y0 = [6 6];
AnonFun = @(t,y) diag([2 + 0.5*y(2) - 1*y(1), 3 + 0.5*y(1) - 1*y(2)])*y;
% AnonFun = @(t,y) [(2 + 0.5*y(2) - 1*y(1))*y(1);
% (3 + 0.5*y(1) - 1*y(2))*y(2)];
[t, y] = ode23(AnonFun, [t0 tfinal], y0);
plot(t, y), grid on
y1 = y(:,1);
y1e = 14/3;
idx = find(t > 1 & y1/y1e > 0.98 & y1/y1e < 1.02); % applying 2% criterion after 1 sec
Ts = t(idx(1))
