I am trying to plot the phase space trajectories for the Hindmarsh-Rose model. I have implemented an RK4 integrator to solve the following set of equations:
The code that I have written so far is given below.
import numpy as np
import matplotlib.pyplot as plt
def RK4(f, x0, t):
dt = t[2] -t[1] #time span
N = len(t)
X = np.empty((len(t), len(x0)))
X[0] = x0
for i in range(1, N):
k1 = f(X[i-1], t[i-1])
k2 = f(X[i-1] + dt/2*k1, t[i-1] + dt/2)
k3 = f(X[i-1] + dt/2*k2, t[i-1] + dt/2)
k4 = f(X[i-1] + dt*k3, t[i-1] + dt)
X[i] = X[i-1] + dt/6*(k1 + 2*k2 + 2*k3 + k4)
return X
def hindmarsh(X, t):
a = 3.0
c = 1.0
d = 5.0
s = 4.0
x0 = - 1.6
# Bifurcation parameters
b = 3.09
I = 3.2
eps = 0.001
x,y,z = X
dxdt = y - (a * x**3) + (b * x**2) + I - z
dydt = c - (d * x**2) - y
dzdt = eps * ( (s * (x - x0)) - z)
return np.array([dxdt, dydt, dzdt])
T = np.linspace(0,100,10000)
Y = [0.03, 0.03, 3]
param = RK4( hindmarsh, Y, T )
ax = plt.axes(projection='3d')
zline = param[2]
yline = param[1]
xline = param[0]
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.plot3D(xline, yline, zline)
However, instead of getting an orbit in the phase space such as the figure below, I get a straight line through the phase space. I would appreciate any tips on how to obtain the plot below.


paramhas shape(len(T), len(Y)), so time is in the first dimension and the x,y,z are in the second dimension. You will get the correct plot with