4次ルンゲクッタ法を実装したいです。

4次ルンゲクッタ法の実装化

pythonで非線形連立微分方程式を4次のルンゲクッタ法で解くことを考えています。
エラーメッセージ等は起こっていないのですが、最終的なグラフの結果が望ましいものとなりません。scipyのodientを用いて解いたものと異なる結果となります。ルンゲクッタ法のプログラムのどこに問題が生じているのか教えて頂けないでしょうか。

発生している問題・エラーメッセージ

scipyのodientを用いた結果と4次ルンゲクッタ法を用いた結果が明らかに異なる。

該当のソースコード

python

4次ルンゲクッタ法を用いたソースコード

import numpy as np
import matplotlib.pyplot as plt
from numpy.typing import _128Bit, _16Bit

パラメータの設定
beta =0.01
lp = 14
ip = 10

区間の分割の設定
T = 100 #
n = 100 # 刻み幅
step = T / n
t = np.arange(0,T,step) #

方程式を定める関数(ラムダ関数)、初期条件の定義
f = lambda S,I,t=0: -betaSI
g = lambda S,E,I,t=0: betaSI-(1/lp)*E
h = lambda E,I,t=0 : (1/lp)*E-(1/ip)*I
q = lambda C,t=0 : (1/ip)*I
S_0 = 1000
E_0 = 10
I_0 = 50
R_0 = 10

結果を返すための配列(行列)の宣言 #初期条件を配列に追加
S = np.empty(n)
E = np.empty(n)
I = np.empty(n)
R = np.empty(n)
S[0] = S_0
E[0] = E_0
I[0] = I_0
R[0] = R_0

方程式を解くための反復計算
def equation():

for i in range(n-1):
    
    _S, _E, _I, _R, _t, = S[i],E[i],I[i],R[i],t[i]
    
    f_1 = step * f(_S,_I,_t)
    f_2 = step * f(_S + f_1/2, _I, _t + step/2)
    f_3 = step * f(_S + f_2/2, _I,_t + step/2)
    f_4 = step * f(_S + f_3, _I, _t + step)
    S[i+1] = _S + 1/6 * (f_1 + 2*f_2 + 2*f_3 + f_4 )

    g_1 = step * g(_S,_E,_I,_t)
    g_2 = step * g(_S + g_1/2, _E, _I, _t + step/2)
    g_3 = step * g(_S + g_2/2, _E, _I, _t + step/2)
    g_4 = step * g(_S + g_3, _E, _I, _t + step )
    E[i+1] = _E + 1/6 * (g_1 + 2*g_2 + 2*g_3 + g_4 )

    h_1 = step * h(_E,_I,_t)
    h_2 = step * h(_E,_I+h_1/2,_t+step/2)
    h_3 = step * h(_E,_I+h_2/2,_t+step/2)
    h_4 = step * h(_E,_I+h_3,_t+step)
    I[i+1] = _I + 1/6 * (h_1 + 2*h_2 + 2*h_3 + h_4 )

    q_1 = step * q(_I,_t)
    q_2 = step * q(_I,_t+step/2)
    q_3 = step * q(_I,_t+step/2)
    q_4 = step * q(_I,_t+step)
    R[i+1] = _R + 1/6 * (q_1 + 2*q_2 + 2*q_3 + q_4 )


return S,E,I,R #

print(I) #データ数はn個
plt.plot(t,S,label = "A",color = "blue")
plt.plot(t,E,label = "B",color = "orange")
plt.plot(t,I,label = "C",color = "green")
plt.plot(t,R,label = "D",color = "red")
plt.legend()
plt.show()


### 試したこと

scipyを用いたソースコード
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as plt

微分方程式の定義
def seir_eq(v,t,beta,lp,ip):
    return [
        -beta*v[0]*v[2],
        beta*v[0]*v[2]-(1/lp)*v[1],
        (1/lp)*v[1]-(1/ip)*v[2],
        (1/ip)*v[2]]

解を求める
ini_state=[1000,10,50,10]
t_max=100
dt=1
t=np.arange(0,t_max,dt)
plt.plot(t,odeint(seir_eq,ini_state,t,args=(0.01,14,10)))

plt.legend(['A','B','C','D'])
plt.show()

こちらの方が望ましい結果が得られます。

### 補足情報（FW/ツールのバージョンなど）

Python 3.8.8