Решение этой системы дифференциальных уравнений с помощью MatLab

Я пытаюсь решить следующую систему дифференциальных уравнений с помощью Matlab.

Параметры, которые я использую, находятся здесь

.

У меня проблема в том, что когда я запускаю свой код, я получаю одинаковые результаты независимо от значений, которые я установил для d1 и d2. Точнее, я получаю те же результаты, что и в упрощенной версии системы, где d1,d2 = 0. Поскольку я пытаюсь повторить чужие результаты, я знаю, что так быть не должно. Кто-нибудь знает, почему это происходит?

% Define u_1(t,j), u_2(t,i) as:
% u_1(t,1) = x(1)
% u_2(t,1) = x(2)
% u_1(t,2) = x(3)
% u_2(t,2) = x(4)

% So

% x(1)' = x(1)*(epsilon - (epsilon/k1)*x(1) - alpha*x(2)) + d1*(rho1(x(4))*x(3) - rho1(x(2))*x(1))
% x(2)' = x(2)*(gama + beta*x(1) - (gama/k2)*x(2)) + d2*(rho2(x(3))*x(4) - rho2(x(1))*x(2))
% x(3)' = x(3)*(epsilon - (epsilon/k1)*x(3) - alpha*x(4)) + d1*(rho1(x(2))*x(1) - rho1(x(4))*x(3))
% x(4)' = x(4)*(gama + beta*x(3) - (gama/k2)*x(4)) + d2*(rho2(x(1))*x(2) - rho2(x(3))*x(4))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

epsilon = 3;
alpha = 0.2;
gama = 0.4;
beta = 0.4;
k1 = 2;
k2 = 1;
m1 = 1;
m2 = 1;
d2 = 0.1;
d1 = 17;

% rho1(x(2)) = m1*exp(-x(2)/m1)
% rho1(x(4)) = m1*exp(-x(4)/m1)
% rho2(x(1)) = m2*exp(-x(1)/m2)
% rho2(x(3)) = m2*exp(-x(3)/m2)

g = @(t,x)[
    x(1)*(epsilon - (epsilon/k1)*x(1) - alpha*x(2)) + d1*(m1*exp(-x(4)/m1)*x(3) - m1*exp(-x(2)/m1)*x(1));
    x(2)*(gama + beta*x(1) - (gama/k2)*x(2)) + d2*(m2*exp(-x(3)/m2)*x(4) - m2*exp(-x(1)/m2)*x(2));
    x(3)*(epsilon - (epsilon/k1)*x(3) - alpha*x(4)) + d1*(m1*exp(-x(2)/m1)*x(1) - m1*exp(-x(4)/m1)*x(3));
    x(4)*(gama + beta*x(3) - (gama/k2)*x(4)) + d2*(m2*exp(-x(1)/m2)*x(2) - m2*exp(-x(3)/m2)*x(4))];

[t1,xa1] = ode45(@(t,x) g(t,x),[0 20],[4 4 4 4]);
[t2,xa2] = ode45(@(t,x) g(t,x),[0 20],[3 3 3 3]);
[t3,xa3] = ode45(@(t,x) g(t,x),[0 20],[2 2 2 2]);
[t4,xa4] = ode45(@(t,x) g(t,x),[0 20],[0.5 0.5 0.5 0.5]);

figure;


subplot(2,2,1)
plot(t1,xa1(:,1),"r",t1,xa1(:,2),"b")
title('y0 = 4')
xlabel('t')
legend({'Free Jobs','Labour Force'},'Location','southwest')

subplot(2,2,2)
plot(t2,xa2(:,1),"r",t2,xa2(:,2),"b")
title('y0 = 3')
xlabel('t')
legend({'Free Jobs','Labour Force'},'Location','southwest')

subplot(2,2,3)
plot(t3,xa3(:,1),"r",t3,xa3(:,2),"b")
title('y0 = 2')
xlabel('t')
legend({'Free Jobs','Labour Force'},'Location','southwest')

subplot(2,2,4)
plot(t4,xa4(:,1),"r",t4,xa4(:,2),"b")
title('y0 = 0.5')
xlabel('t')
legend({'Free Jobs','Labour Force'},'Location','southwest')
sgtitle('Labour Force and Free Jobs')


figure;

hold on

plot(t1,xa1(:,1),"r",t1,xa1(:,2),"b")
plot(t2,xa2(:,1),"r",t2,xa2(:,2),"b")
plot(t3,xa3(:,1),"r",t3,xa3(:,2),"b")
plot(t4,xa4(:,1),"r",t4,xa4(:,2),"b")
title('Labour Force and Free Jobs')
xlabel('t')
legend({'Free Jobs','Labour Force'},'Location','northeast')

hold off