title: Solving this system of differential equations with MatLab
body: <p>I'm trying to solve the following system of differential equations with Matlab.</p>
  
  <p><img src="https://i.stack.imgur.com/6Yncj.png" alt="the following system"></p>
  
  <p>The parameters i use are found here</p>
  
  <p><img src="https://i.stack.imgur.com/6e88E.png" alt="are found here">. </p>
  
  <p>The problem i have is that when i run my code i get the same results regardless of the values i set for d1 and d2. Precisely, i get the same results as in the simplified version of the system, where d1,d2 = 0. Since i'm trying to replicate the results of someone else, i know that this should not be the case. Does anyone have any idea why this is happening?</p>
  
  <pre><code>% 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
  </code></pre>