在k1、k2和k3的计算中，您有一个额外的f(x0)调用。将函数kcontants更改为

def kcontants(f,h,x0):
    k0=h*f(x0)
    k1=h*f(x0 + k0/2)
    k2=h*f(x0 + k1/2)
    k3=h*f(x0 + k2)
    #note returned K is a matrix
    return np.array([k0,k1,k2,k3])

你是否考虑过不同的初始值来进行计算？你选择的初始值是否合理？也就是说，它们是物理的吗？从以往的经验来看，如果你选择了愚蠢的初始参数，你有时会得到非常混乱的结果。

晚安。这是我用scipyedo集成器scipy.integrate.odeint做的版本。

# Author      : Carlos Eduardo da Silva Lima
# Theme       : Movement of a Plant around a fixed star
# Language    : Python
# date        : 11/19/2022
# Environment : Google Colab

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import root
from scipy.linalg import eig
from mpl_toolkits.mplot3d import Axes3D

##################################
# Condições inicial e parãmetros #
##################################
t_inicial = 0
t_final = 100
N = 10000
h = 1e-3

x_0 = 1.0
y_0 = 1.0
z_0 = 1.0

#####################
# Equação de Lorenz #
#####################
def Lorenz(r,t,sigma,rho,beta):

  x = r[0]; y = r[1]; z = r[2]

  edo1 = sigma*(y-x)
  edo2 = x*(rho-z)-y
  edo3 = x*y-beta*z
  return np.array([edo1,edo2,edo3])

t = np.linspace(t_inicial,t_final,N)
r_0 = np.array([x_0,y_0,z_0])
#sol = odeint(Lorenz,r_0,t,rtol=1e-6,args = (10,28,8/3))
sol = odeint(Lorenz, r_0, t, args=(10,28,8/3), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=1e-9, atol=1e-9, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False)

'''x = sol[:,0]
y = sol[:,1]
z = sol[:,2]'''
x, y, z = sol.T

# Plot
plt.style.use('dark_background')
ax = plt.figure(figsize = (10,10)).add_subplot(projection='3d')
ax.plot(x,y,z,'m-',lw=0.5, linewidth = 1.5)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title("Atrator de Lorenz")
plt.show()

在第二部分中，我模拟了两个Lorenz系统，以验证系统对初始条件的敏感依赖性，在第二个系统中，我在x（t0），y（t0）和z（t0）的初始条件上加上了一定量的eps = 1 e-3。

# Depedência com as condições iniciais
eps = 1e-3
r_0_eps = np.array([x_0+eps,y_0+eps,z_0+eps])
sol_eps = odeint(Lorenz, r_0_eps, t, args=(10,28,8/3), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None,
             rtol=1e-9, atol=1e-9, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False)

'''x_eps = sol_eps[:,0]
y_eps = sol_eps[:,1]
z_eps = sol_eps[:,2]'''
x_eps, y_eps, z_eps = sol_eps.T

# Plot
plt.style.use('dark_background')
ax = plt.figure(figsize = (10,10)).add_subplot(projection='3d')
ax.plot(x,y,z,'r-',lw=1.5)
ax.plot(x_eps,y_eps,z_eps,'b-.',lw=1.1)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title("Lorenz Attractor")
plt.show()

希望我能帮上忙，再见：）。