SageMathでは，GSLを使って微分方程式を解くこともできる．desolve_system_rk4よりも速いようだ．まずは$\sin$と$\cos$が解となる微分方程式 \[ \begin{aligned} \frac{dy_0}{dt} & = y_1,\\ \frac{dy_1}{dt} & = -y_0 \end{aligned} \] を$0\le t \le 10$の範囲で解いてみる．

T = ode_solver() # 微分方程式を解くオブジェクトを得る
T.algorithm = "rk4" # 四次のRunge–Kutta法を使う
# 微分方程式系の右辺を返す関数を定義．yは関数の値で，未知関数の分だけのリストで与えられる．リストなので0始まり．
def f(t,y):
    return [y[1],-y[0]]

T.function = f # 微分方程式系の右辺をセット
T.y_0 = [0,1] # 未知関数の初期値．それぞれy_0(0),y_1(0)．tの初期値は後で与える．
# 微分方程式を解く．t_spanが変数tの範囲で，num_points=100はそこを100分割して計算する．
T.ode_solve(t_span = [0,10],num_points = 100)
# 解はT.solutionに入っている．(t_n,[y_0(t_n),y_1(t_n)])というものからなる配列．desolve_system_rk4とちょっと違う．
list_plot([(t_1,y_1[0]) for t_1,y_1 in T.solution], plotjoined=True)
# こういうのもある．何かグラフがでる．
T.plot_solution()

t_spanとnum_pointsを与える代わりに，t_spanに点の配列を与えるのでもよいらしい．上は次と（ほぼ）等価．

T.ode_solve(t_span = [x/10.0 for x in range(101)])

アルゴリズムによっては微分方程式の右辺を未知関数で偏微分したものが必要になることがある．微分方程式系が \[ \begin{aligned} \frac{dy_0}{dt} & = f_0(t,y_0,\ldots,y_n),\\ \frac{dy_1}{dt} & = f_2(t,y_1,\ldots,y_n),\\ \cdots\\ \frac{dy_n}{dt} & = f_n(t,y_1,\ldots,y_n) \end{aligned} \] で与えられているならば，T.jacobianに \[ [ [\frac{\partial f_0}{\partial y_0},\frac{\partial f_0}{\partial y_1},\ldots,\frac{\partial f_0}{\partial y_n}], ..., [\frac{\partial f_n}{\partial y_0},\frac{\partial f_n}{\partial y_1},\ldots,\frac{\partial f_n}{\partial y_n}], [\frac{\partial f_0}{\partial t},\frac{\partial f_1}{\partial t},\ldots,\frac{\partial f_n}{\partial t}] ] \] という値を返す関数を設定しておく．（関数は上のT.functionと同様にPythonの関数として定義する．）

T.functionとT.jacobianに設定する関数はさらにもう一つ引数をとることができる．これはパラメータを表す．前のジャパニーズ・アトラクタの例 \[ \begin{aligned} \frac{dy_0}{dt} & = y_1,\\ \frac{dy_1}{dt} & = -ky_1 - y_0^3 + B\cos(t) \end{aligned} \] という方程式（$k$，$B$がパラメータ）の場合は次のようになる．

def f(t,y,params):
  return [y[1],-params[0]*y[1] - y[0]^3 + params[1]*cos(t)]

T.function = f
def j(t,y,params):
  return [[0,1],[-3*y[0]^2,-params[0]],[0,-params[1]*sin(t)]]

T.jacobian = j
# デフォルトのアルゴリズムで（なのでjacobianは不要だけど……）
T.algorithm = "rkf45"
# params=でパラメータを与える．
T.ode_solve(t_span = [pi*x/30 for x in range(300001)],params=[0.05,7])
list_plot([(T.solution[60*i][1][0],T.solution[60*i][1][1]) for i in range(len(T.solution)/60)])

時間もとってみたよ．

CPU times: user 4min 27s, sys: 28.9 s, total: 4min 55s
Wall time: 4min 58s

話にならないなdesolve_system_rk4．