python　ルンゲクッタ法

python
1import numpy as np
2import matplotlib.pyplot as plt
3from numpy.lib.type_check import isrealobj
4
5#関数の定義
6def f1(fai2):
7    return fai2
8
9def f2(fai1,fai2):
10    return (i-np.sin(fai1)-fai2)/bc
11
12#初期値
13a = 0.0
14b = 2.0
15N = 10000
16h = (b-a)/N
17j = 0
18
19i = 1.2
20bc = 2.0
21
22t = np.arange(a,b,h)
23fai1_0 = 0.0
24fai2_0 = 0.0
25
26#配列
27fai1 = np.empty(N)
28fai2 = np.empty(N)
29fai1[0] = fai1_0
30fai2[0] = fai2_0
31
32for j in range (N-2):
33    k1 = h * f1(fai2[j])
34    d1 = h * f2(fai1[j] , fai2[j])
35    k2 = h * f1(fai2[j] + d1/2)
36    d2 = h * f2(fai1[j] + k1/2, fai2[j] + d1/2)
37    k3 = h * f1(fai2[j] + d2/2)
38    d3 = h * f2(fai1[j] + k2/2, fai2[j] +d2/2)
39    k4 = h * f1(fai2[j] + d3)
40    d4 = h * f2(fai1[j] + k3 , fai2[j] + d3)     
41    fai1[j+1] = fai1[j] + 1/6 * (k1 + 2 * k2 + 2 * k3 + k4) 
42    fai2[j+1] = fai2[j] + 1/6 * (d1 + 2 * d2 + 2 * d3 + d4) 
43
44    if j == N-1:
45        break
46
47    j += 1
48
49
50v = np.average(fai2[6000:])
51print(v)
52print(i)
53plt.scatter (v, i, marker='.') 
54
55#plot
56plt.plot(v,i)   
57plt.xlim(0, 1.8)
58plt.ylim(0, 2.0)
59plt.xlabel("v")
60plt.ylabel("i")
61plt.show()
62```ルンゲクッタ法を用いたコードです。このコードを実行すると、![![イメージ説明](3252b50fc1815f3d04ce8dfc98d91da6.png)](f7cee8d6b034ad803dd0ce8d10181ea1.png)
63となるのですが、本来であれば、
64![イメージ説明](9418446d9abaf6a6fcec677c60b2a796.jpeg)
65のように、i=1.2,Bc=2.0のときは、この点はv=1.2～1.4付近にプロットされなければなりません。
66どこがおかしいでしょうか。