Hi Mohamed H Hassan

I read this Paper and I Hope this code can Help you :

% Define the DDEs

function dydt = sir_model(t, y, Z)

% Parameters from paper

beta = ...; % Fill in from paper

gamma = ...; % Fill in from paper

mu = ...; % Fill in from paper

N = ...; % Fill in from paper

% ... Add other parameters

% Previous values

S_tau = Z(1,:);

I_tau = Z(2,:);

% Equations from paper

dS = ...; % Fill in from paper

dI = ...; % Fill in from paper

dR = ...; % Fill in from paper

dydt = [dS; dI; dR];

end

% Define delays

lags = ...; % Fill in from paper

% Initial conditions and history

y0 = ...; % Fill in from paper

history = ...; % Fill in from paper

% Time span

tspan = [0 100]; % Example

% Solve using dde23

sol = dde23(@sir_model, lags, y0, tspan, [], history);

% Plot the results

figure;

plot(sol.x, sol.y(1,:), '-r', sol.x, sol.y(2,:), '-g', sol.x, sol.y(3,:), '-b');

legend('Susceptible', 'Infected', 'Recovered');

xlabel('Time');

ylabel('Population');

title('Delayed SIR Model');

mathworks.com

see this

https://www.mathworks.com/matlabcentral/answers/460308-how-to-plot-system-of-fractional-delay-differential-equations