def rk(func, t0, init, h): f1 = func(t0 , init ) f2 = func(t0 + h / 2, init + h * f1 / 2) f3 = func(t0 + h / 2, init + h * f2 / 2) f4 = func(t0 + h , init + h * f3 ) return init + (h / 6) * (f1 + 2 * (f2 + f3) + f4) # 次の微分方程式を 1/8 刻みで 1 ステップだけ解いてみる。真値は 0.236111111...... # x'(t) = x(t)^2 - t^2 - 2t + 2 # x(0) = 0 def xDot(t, x): return x**2 - t**2 - 2 * t + 2 print(rk(xDot, 0, 0, 1/8))