I presume that Xinit is not equal to 
, and 
. I tested this with the ode45 solver, and it also works with non-zero initial values.
, and . I tested this with the ode45 solver, and it also works with non-zero initial values.% Parameters
Xinit = 1;
frac = 1/2;
rext = 1;
rtra = 1;
rvol = 1;
Sr = 1;
Sti = 1;
Sfeu = 1;
% Create a function handle F for a system of 1st-order ODEs:
F = @(t, y) [- 0.000038*y(1) - (y(1)*(y(1)/Xinit)^frac)*rext;
- 0.000038*y(2) + rext*y(1) - rtra*y(2) + Sr;
- 0.000038*y(3) + rext*y(2) - rtra*y(3) + Sti;
0.000038*y(4) + rext*y(3) - rvol*y(4) + Sfeu];
tspan = [0 10];
y0 = [3; 2; 1; 0];
[t, y] = ode45(F, tspan, y0);
% Plotting the result
plot(t, y, "-o"), grid on
xlabel('Time, t (seconds)'), ylabel('\bf{y}(t)')
legend('X', 'Y', 'Z', 'U', 'location', 'NW')
% Parameters
Xinit = 1;
frac = 1/2;
rext = 1;
rtra = 1;
rvol = 1;
Sr = 1;
Sti = 1;
Sfeu = 1;
% Setting up the ODE object:
F = ode;
F.InitialValue = [3; 2; 1; 0];
F.ODEFcn = @(t, y) [- 0.000038*y(1) - (y(1)*(y(1)/Xinit)^frac)*rext;
- 0.000038*y(2) + rext*y(1) - rtra*y(2) + Sr;
- 0.000038*y(3) + rext*y(2) - rtra*y(3) + Sti;
0.000038*y(4) + rext*y(3) - rvol*y(4) + Sfeu];
F.Solver = "ode45";
S = solve(F, 0, 10); % Solving F over the time from 0 to 10 s
% Plotting the result
plot(S.Time, S.Solution, "-o"), grid on
xlabel('Time, t (seconds)'), ylabel('\bf{y}(t)')
legend('X', 'Y', 'Z', 'U', 'location', 'NW')
