Matlab coding ordinary differential equation initial condition graph not running in the initial condition

Matlab coding ordinary differential equation initial condition graph not running in the initial condition. For x(2), x(6) and x(7) only getting properly graph is coming with the initial condition starting but for x(1), x(3), x(4), x(5) graph is not coming as the initial condition starting

Tr()

30 successful steps
2 failed attempts
193 function evaluations
Elapsed time is 0.072265 seconds.

function Tr
options = odeset('RelTol',1e-6,'Stats','on');
%initial conditions    
Xo = [0.005; 0.0007; 0.001; 0.0001; 0.0001; 0; 0];  
tspan = [0,120];
tic
[t,X] = ode45(@TestFunction,tspan,Xo,options);
toc  
%figure
hold on
plot(t, X(:,1), 'green')
plot(t, X(:,2), 'blue')
plot(t, X(:,3), 'red')
plot(t, X(:,4), 'yellow')
plot(t, X(:,5), 'green')
plot(t, X(:,6), 'cyan')
plot(t, X(:,7), 'y')
%hold on 
% legend('x1','x2')
% ylabel('x - Population')
% xlabel('t - Time')
%hold on
end
function [dx_dt]= TestFunction(~,x)
% Parameters
s       = 0.038;
alpha   = 0.02;
gamma   = 0.10;
r       = 0.03;
dH      = 0.0083;
dX      = 0.0125;
rho     = 0.07;
K       = 500;
beta    = 0.0005;
theta   = 0.03;
eta     = 0.015;
muH     = 0.015;
muI     = 0.08;
lambdaP = 0.05;
lambdaL = 0.043;
UP      = 0.20;
deltaP  = 0.033;
UL      = 0.50;
deltaL  = 0.05;
kappa   = 0.005;
dx_dt(1) = (s) - ((alpha * x(3) * x(1)) / (1 + gamma * x(1))) - (dH * x(1)) + (r * x(2));   
dx_dt(2) = ((alpha * x(3) * x(1)) / (1 + gamma * x(1))) - ((r+dX) * x(2));   
dx_dt(3) = ((rho * x(3)) * (1 - (x(3) + x(4)) / K)) - (beta * x(3) * x(5)) - (lambdaP * x(6) * x(3)) - (lambdaL * x(7) * x(3)) - (muH * x(3)); 
dx_dt(4) = (beta * x(3) * x(5)) - (lambdaP * x(6) * x(4)) - (lambdaL * x(7) * x(4)) - (muI * x(4)) - (kappa * x(5) * x(4));
dx_dt(5) = (theta * x(4)) - (eta * x(5));
dx_dt(6) = (UP) - (deltaP * x(6));      
dx_dt(7) = (UL) - (deltaL * x(7));
dx_dt = dx_dt';  
end
%hold on

[Torsten]

If the results for x(1), x(3), x(4) and x(5) make no sense for you, you will have to check your parameters and equations for dx(1), dx(3), dx(4) and dx(5).

[Sam Chak]

Hi @Dhivyadharshini

The responses for x1, x3, x4, and x5 begin from the specified initial values: {0.005, 0.001, 0.0001, 0.0001}.

options = odeset('RelTol', 1e-6);
X0      = [0.005; 0.0007; 0.001; 0.0001; 0.0001; 0;  0];  
        %  x1     x2      x3     x4      x5      x6  x7
tspan   = [0, 120];
[t, X]  = ode45(@TestFunction, tspan, X0, options);

% figure(1)
% tL  = tiledlayout(2, 2, 'TileSpacing', 'Compact');
% 
% nexttile
% plot(t, X(:,1), 'color', '#AA0815'), grid on 
% title('x1')
% 
% nexttile
% plot(t, X(:,3), 'color', '#F0B41C'), grid on 
% title('x3')
% 
% nexttile
% plot(t, X(:,4), 'color', '#49B6A9'), grid on 
% title('x4')
% 
% nexttile
% plot(t, X(:,5), 'color', '#3D9BE1'), grid on 
% title('x5')
% 
% xlabel(tL, 'Time')
% ylabel(tL, 'Population')

figure(2)
tL2 = tiledlayout(2, 2, 'TileSpacing', 'Compact');

nexttile
plot(t, X(:,1), 'color', '#AA0815'), grid on 
xlim([0, 120])
title('x1')

nexttile
plot(t, X(:,3), 'color', '#F0B41C'), grid on 
xlim([0, 120])
title('x3')

nexttile
plot(t, X(:,4), 'color', '#49B6A9'), grid on 
xlim([0, 120])
title('x4')

nexttile
plot(t, X(:,5), 'color', '#3D9BE1'), grid on 
xlim([0, 120])
title('x5')

%title(tL2, 'Zoomed-in')
xlabel(tL2, 'Time')
ylabel(tL2, 'Population')

%% Dynamics of Population Growth
function [dx_dt]= TestFunction(~,x)
    % Parameters
    s       = 0.038;
    alpha   = 0.02;
    gamma   = 0.10;
    r       = 0.03;
    dH      = 0.0083;
    dX      = 0.0125;
    rho     = 0.07;
    K       = 500;
    beta    = 0.0005;
    theta   = 0.03;
    eta     = 0.015;
    muH     = 0.015;
    muI     = 0.08;
    lambdaP = 0.05;
    lambdaL = 0.043;
    UP      = 0.20;
    deltaP  = 0.033;
    UL      = 0.50;
    deltaL  = 0.05;
    kappa   = 0.005;
    dx_dt(1) = (s) - ((alpha * x(3) * x(1)) / (1 + gamma * x(1))) - (dH * x(1)) + (r * x(2));   
    dx_dt(2) = ((alpha * x(3) * x(1)) / (1 + gamma * x(1))) - ((r+dX) * x(2));   
    dx_dt(3) = ((rho * x(3)) * (1 - (x(3) + x(4)) / K)) - (beta * x(3) * x(5)) - (lambdaP * x(6) * x(3)) - (lambdaL * x(7) * x(3)) - (muH * x(3)); 
    dx_dt(4) = (beta * x(3) * x(5)) - (lambdaP * x(6) * x(4)) - (lambdaL * x(7) * x(4)) - (muI * x(4)) - (kappa * x(5) * x(4));
    dx_dt(5) = (theta * x(4)) - (eta * x(5));
    dx_dt(6) = (UP) - (deltaP * x(6));      
    dx_dt(7) = (UL) - (deltaL * x(7));
    dx_dt = dx_dt';  
end

[Dhivyadharshini]

I want for 120 days x-axis 120 should come

[Torsten]

Then adapt the "xlim" commands in the code (see above).