Finite difference method converging to wrong steady state (python)

I'm trying to recreate numerical results from an article. In particular, I want to numerically solve the following PDE system:

In this system, Fn represents a prey population density, Cn a predator population density, and En the energy of the predators. The (3.57) equations model how these different populations interact when time t lies between two consecutive natural numbers, i.e., n<t<n+1 for some natural number n. The (3.4) equations govern how the populations change when time is a natural number, i.e., when t=n for some natural number n.

When the authors numerically solve the non-spatial system, i.e., with all the diffusion terms in the (3.57) equations dropped, they produce the following results (the parameters they use are at the bottom of the image):

When I numerically solve the system in Python, I get the following results:

My results look very similar, however, the peaks of my solutions are always around twice as high as the solutions from the article. This discrepancy doesn't disappear by taking smaller time steps so I'm thinking something is wrong with how I'm implementing the forward Euler method in Python. Any help would be greatly appreciated.

import matplotlib.pyplot as plt

ρ=0.55
r1=10
θ1=0.17
θ2=0.1
m1=2.1
β=7
ω=0.01
δ1=0.05

h=0.01 #step size
T=15000 #number of steps

D=[] #domain

for i in range(0,T+1):
    D.append(i*h)

F=0.05 #prey initial condition
C=2.5 #predator initial condition
E=0 #energy initial condition

Fval=[F] #list that will store all calculated prey values
Cval=[C] #list that will store all calculated predator values
Eval=[E] #list that will store all calculated energy values

for i in range(1,T+1):
    if (i*h)%1==0: #between year dynamics
        Fnew=ρ*F
        Cnew=E+C
        Enew=0
    else: #within year dynamics
        Fnew=(1+h*r1*(1-F))*F-((h*F*C)/(ω+θ1*F+θ2*C))
        Cnew=(1-h*m1)*C
        Enew=((h*β*F*C)/(ω+θ1*F+θ2*C))+(1-h*δ1)*E
        
    Fval.append(Fnew)
    Cval.append(Cnew)
    Eval.append(Enew)

    F=Fnew
    C=Cnew
    E=Enew

#plotting
plt.plot(D, Fval,color='red',linestyle='dotted')
plt.plot(D, Cval,color='blue',linestyle='dashdot')
plt.plot(D, Eval,color='black')
plt.xlabel('time t')
plt.ylabel('$F_n(t)$, $C_n(t)$, $E_n(t)$')
plt.legend(['$F_n(t)$', '$C_n(t)$','$E_n(t)$']) 
plt.xlim(0,50)
plt.ylim(0)
plt.tick_params(axis="x", direction="in")
plt.tick_params(axis="y", direction="in")
plt.title('Spatial homogeneous impulsive solutions')
plt.show()

plt.plot(D, Fval,color='blue')
plt.xlabel('time t')
plt.ylabel('$F_n(t)$')
plt.xlim(140,150)
plt.ylim(0)
plt.tick_params(axis="x", direction="in")
plt.tick_params(axis="y", direction="in")
plt.title('Spatial impulsive solution')
plt.show()

plt.plot(D, Cval,color='blue')
plt.xlabel('time t')
plt.ylabel('$C_n(t)$')
plt.xlim(140,150)
plt.ylim(0)
plt.tick_params(axis="x", direction="in")
plt.tick_params(axis="y", direction="in")
plt.title('Spatial impulsive solution')
plt.show()

plt.plot(D, Eval,color='blue')
plt.xlabel('time t')
plt.ylabel('$E_n(t)$')
plt.xlim(140,150)
plt.ylim(0)
plt.tick_params(axis="x", direction="in")
plt.tick_params(axis="y", direction="in")
plt.title('Spatial impulsive solution')
plt.show()

Well, I tried coding it up independently. Sorry, it's in C++, with graphs from gnuplot. It's also using 4th-order Runge-Kutta.

Basically, ... I agree with YOUR results!

The final state doesn't seem sensitive to initial conditions, so I think the authors were perhaps using different values of the coefficients - perhaps you could try changing them.

Since at least this gives you independent confirmation of your result it might be worth contacting the authors and asking them to check their coefficients.

#include <iostream>
#include <valarray>
using namespace std;

const double rho = 0.55, r1 = 10.0, theta1 = 0.17, theta2 = 0.1, m1 = 2.1, b = 7, omega = 0.01, delta1 = 0.05;

//======================================================================

valarray<double> rk4( double x, const valarray<double> &Y, valarray<double> f( double, const valarray<double> & ), double dx )
{
   valarray<double> dY1, dY2, dY3, dY4;
   dY1 = dx * f( x           , Y             );
   dY2 = dx * f( x + 0.5 * dx, Y + 0.5 * dY1 );
   dY3 = dx * f( x + 0.5 * dx, Y + 0.5 * dY2 );
   dY4 = dx * f( x +       dx, Y +       dY3 );
   return Y + ( dY1 + 2.0 * dY2 + 2.0 * dY3 + dY4 ) / 6.0;
}

//======================================================================

valarray<double> deriv( double t, const valarray<double> &Y )
{
   double F = Y[0], C = Y[1], E = Y[2];
   double dFdt = r1 * ( 1.0 - F ) * F - F * C / ( omega + theta1 * F + theta2 * C );
   double dCdt = -m1 * C;
   double dEdt = b * F * C / ( omega + theta1 * F + theta2 * C ) - delta1 * E;
   return valarray<double>{ dFdt, dCdt, dEdt };
}

//======================================================================

int main()
{
   valarray<double> Y = { 0.05, 2.5, 0.0 };
   const int N1 = 100;   // timesteps per year
   const double dt = 1.0 / N1;
   const int Years = 150.0;
   double t = 0;
   for ( int y = 0; y < Years; y++ )
   {
      // Within a year
      for ( int n = 0; n < N1; n++ )
      {
         cout << t << '\t' << Y[0] << '\t' << Y[1] << '\t' << Y[2] << '\n';
         Y = rk4( t, Y, deriv, dt );
         t += dt;
      }
      // End-of-year reset
      Y[0] = rho * Y[0];
      Y[1] = Y[2] + Y[1];
      Y[2] = 0.0;
   }
}

Graphs