Finite Difference Method implementation in Python

Question

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. Note zero-flux boundary conditions are imposed, i.e.,

and similarly for Cn and En.

When the authors numerically solve the system (they claim to use a finite difference scheme), with the following parameters,

they produce the following results:

When I run my simulation with D1=10, however, I get the following:

Which looks like what the others got for D1=5. Is there something wrong with my implementation of the finite difference method? Am I implementing the zero-flux boundary conditions wrong? Any help would be greatly appreciated.

from matplotlib.animation import FuncAnimation
import numpy as np
import matplotlib.pyplot as plt
import random

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

k=0.01 #time step
K=5000 #number of time steps

h=0.8 #space step
L=20 #length of space domain
H=int(L/h) #number of space steps

print('r =',(D1*k)/(h**2))

X=np.linspace(0,L,H+1) #space domain

T=np.linspace(0,K*k,K+1) #time domain


F=[] #prey initial conditions
for i in range(0,H+1):
    F.append(0.03+(0.01*(random.random())))

C=[] #predator initial conditions
for i in range(0,H+1):
    C.append(2.4+(0.05*(random.random())))

E=[] #energy initial conditions
for i in range(0,H+1):
    E.append(0.02+(0.05*(random.random())))


Fval=[] #prey values
for i in range(0,K+1):
    Fval.append([])

Cval=[] #predator values
for i in range(0,K+1):
    Cval.append([])

Eval=[] #energy values
for i in range(0,K+1):
    Eval.append([])

Fval[0]=F
Cval[0]=C
Eval[0]=E

for i in range(1,K+1):
    if (i*k)%1==0: #between year dynamics
        Fval[i]=[x*ρ for x in Fval[i-1]]
        Cval[i]=[x+y for x,y in zip(Eval[i-1],Cval[i-1])]
        Eval[i]=[x*0 for x in Eval[i-1]]
        F=Fval[i]
        C=Cval[i]
        E=Eval[i]
        continue
    for j in range(0,H+1): #within year dynamics
        if j==0: #zero flux boundary condition 
            Fnew=(1+k*r1*(1-F[j]))*F[j]-((k*F[j]*C[j])/(ω+θ1*F[j]+θ2*C[j]))+((2*k)/(h**2))*(F[j+1]-F[j])
            Cnew=(1-k*m1)*C[j]+D1*((2*k)/(h**2))*(C[j+1]-C[j])
            Enew=((k*β*F[j]*C[j])/(ω+θ1*F[j]+θ2*C[j]))+(1-k*δ1)*E[j]+D1*((2*k)/(h**2))*(E[j+1]-E[j])
        elif j==H: #zero flux boundary condition
            Fnew=(1+k*r1*(1-F[j]))*F[j]-((k*F[j]*C[j])/(ω+θ1*F[j]+θ2*C[j]))+((2*k)/(h**2))*(F[j-1]-F[j])
            Cnew=(1-k*m1)*C[j]+D1*((2*k)/(h**2))*(C[j-1]-C[j])
            Enew=((k*β*F[j]*C[j])/(ω+θ1*F[j]+θ2*C[j]))+(1-k*δ1)*E[j]+D1*((2*k)/(h**2))*(E[j-1]-E[j])
        else:
            Fnew=(1+k*r1*(1-F[j]))*F[j]-((k*F[j]*C[j])/(ω+θ1*F[j]+θ2*C[j]))+(k/(h**2))*(F[j+1]-2*F[j]+F[j-1])
            Cnew=(1-k*m1)*C[j]+D1*(k/(h**2))*(C[j+1]-2*C[j]+C[j-1])
            Enew=((k*β*F[j]*C[j])/(ω+θ1*F[j]+θ2*C[j]))+(1-k*δ1)*E[j]+D1*(k/(h**2))*(E[j+1]-2*E[j]+E[j-1])
        Fval[i].append(Fnew)
        Cval[i].append(Cnew)
        Eval[i].append(Enew)
    F=Fval[i]
    C=Cval[i]
    E=Eval[i]

#plotting
norm = plt.Normalize(np.min(Cval), np.max(Cval))

ax = plt.axes(projection='3d')
for i in range(0,K+1):
    ax.scatter3D(X,[i*k],Cval[i],c=Cval[i],cmap='inferno',norm=norm);
ax.set_xlabel('space x')
ax.set_ylabel('time t')
ax.set_zlabel('$C_n(x,t)$')
ax.xaxis.set_rotate_label(False)
ax.yaxis.set_rotate_label(False)
ax.zaxis.set_rotate_label(False)
ax.view_init(azim=-120)
plt.show()