Looks a lot like you have (at least one) bifurcation point.

You can perform a stability analysis based on linerization which shows that actually both equilibria have exactly one unstable direction corresponding to an eigenvalue with negative real part.

Thus, if the trajectories from the different solvers are somewhat sufficiently aligned with those unstable directions, this causes them to miss that equilibrium / steady state. Conversely, if they progress in directions more or less aligned with the stable ones, they approach the corresponding equilibrium.
In the context of your system, `ode23s`

"finds" thus the stable directions of the $(1, 0, 0, 0, 0)$ equilibrium or is repelled by the unstable direction of the $\sim (0.0589, 0.0021, 0.7937, 0.0655, 0.0799)$ steady state (or mixture of both).
Converse holds for the `ode23tb`

solver.

I attach some Matlab files I used throughout the analysis, including a animation of the trajectories in reduced ($I, L, R + L +W$) space, showing that somewhere around $I \approx L \approx 0, R + L + W \approx 1$ the bifurcation happends.

```
function f = ODE_System(~, SIRLW)
% Parameters
alpha = 26;
beta = 260;
kappa = 0.1;
k = 3 * kappa; % Kappa only used after multiplication with 3
gamma = 17;
nu = 5;
xi = 0.0125;
S = SIRLW(1); I = SIRLW(2); R = SIRLW(3); L = SIRLW(4); W = SIRLW(5);
f = zeros(5, 1);
f(1) = -beta * I * S + xi * (1 - S) + k * L;
f(2) = beta * I * S + alpha * I * L - xi * I - gamma * I;
f(3) = nu * beta * I * W + gamma * I - xi * R - k * R;
f(4) = -alpha * I * L + k * W - xi * L - k * L;
f(5) = - nu * beta * I * W + k * R - xi * W - k * W;
end
```

```
function Df_Dx = Eval_Jacob(SIRLW_eq)
% Evaluates Jacobian at given point
% Parameters
alpha = 26;
beta = 260;
kappa = 0.1;
k = 3 * kappa; % Kappa only used after multiplication with 3
gamma = 17;
nu = 5;
xi = 0.0125;
S_eq = SIRLW_eq(1); I_eq = SIRLW_eq(2); R_eq = SIRLW_eq(3); L_eq = SIRLW_eq(4); W_eq = SIRLW_eq(5);
Df_Dx = [-beta * I_eq - xi, -beta * S_eq, 0, k, 0;
beta * I_eq, beta * S_eq + alpha * L_eq - xi - gamma, 0, alpha * I_eq, 0;
0, nu * beta * W_eq + gamma, -xi - k, 0, nu * beta * I_eq;
0, -alpha * L_eq, 0, -alpha * I_eq - xi - k, 0;
0, -nu * beta * W_eq, k, 0, -xi - k; ];
end
```

```
clear;
close all;
clc;
% Parameters
alpha = 26;
beta = 260;
kappa = 0.1;
k = 3 * kappa; % Kappa only used after multiplication with 3
gamma = 17;
nu = 5;
xi = 0.0125;
%% Equilibrium of upper / first plot
SILRW_Eq_1 = [1; 0; 0; 0; 0];
Df_Dx = Eval_Jacob(SILRW_Eq_1);
[EVec_Eq_1, EVal_Eq_1] = eigs(Df_Dx); % Shows one eigenvalue (whose real part) > 0 => Unstable equilibrium
%% Equilibrium of lower / second plot
tspan = [0, 100];
S_0 = 0.99;
SILRW_0 = [S_0; 1 - S_0; 0; 0; 0];
[t_tb, SILRW_tb] = ode23tb(@ODE_System, tspan, SILRW_0);
[t_s, SILRW_s] = ode23s( @ODE_System, tspan, SILRW_0);
% Show that this values lead indeed to (more or less) an equilibrium:
SILRW_Eq_2 = SILRW_tb(end, :)';
t_dummy = 42;
disp(ODE_System(t_dummy, SILRW_Eq_2)); % Not really 0 yet
% Other way of obtaining this equilibrium: Solve eqs. numerically with ode solver outcome as start value
options = optimoptions('fsolve', 'Display', 'iter');
fun = @(x) ODE_System(t_dummy, x);
SILRW_Eq_2 = fsolve(fun, SILRW_Eq_2, options);
% Check values of RHS again, they are now machine precision:
disp(ODE_System(t_dummy, SILRW_Eq_2));
Df_Dx = Eval_Jacob(SILRW_Eq_2);
[EVec_Eq_2, EVal_Eq_2] = eigs(Df_Dx); % Shows one eigenvalue with real part > 0 => Unstable equilibrium
%% Investigate trajectories
% Bundle RLW into one variable (leaves one with 3 variables which can be plotted)
SI_Bundled_tb = [SILRW_tb(:, 1), SILRW_tb(:, 2), SILRW_tb(:, 3) + SILRW_tb(:, 4) + SILRW_tb(:, 5)];
SI_Bundled_s = [SILRW_s(:, 1), SILRW_s(:, 2), SILRW_s(:, 3) + SILRW_s(:, 4) + SILRW_s(:, 5)];
for i = 1:length(t_s)
plot3(SI_Bundled_s(i, 1), SI_Bundled_s(i, 2), SI_Bundled_s(i, 3), '-rx');
xlabel('S');
ylabel('I');
zlabel('R + L + W');
title('Trajectories. Red: ode23 S, Blue: ode23 TB');
hold all;
pause(0.1);
end
for i = 1:length(t_tb) / 3 % All the action done by then
plot3(SI_Bundled_tb(i, 1), SI_Bundled_tb(i, 2), SI_Bundled_tb(i, 3), '-bo', 'Displayname', 'ode23_TB');
hold all;
pause(0.1);
end
```