関数を変えるとコードの結果が正しく得られない

実現したいこと

結果が正しく出るようにコードを直したい

前提

pythonで３つの接合が並列接合された際の数値解析を調べています。
印加する電流を直流から交流に変えると、本来得られるはずのグラフとは異なった結果が出てきます。
直流の時は正しい結果が出ていたので、おかしいのは交流に変更した際の関数の表記の仕方だと思い、時間の表記を変えたり、コードを細かく書いたりして試行錯誤を繰り返したのですが、結果はずっと同じままでした。
ぜひ改善点をお教えいただきたいです。

該当のソースコード

import numpy as np import matplotlib.pyplot as plt import math pi = math.pi faia11 = float() faia21 = float() faia12 = float() faia22 = float() faib1 = float() faib2 = float() tau = float() I = float() t = float() A = 1 w = 0.01 def f1_1(tau,faia21): return faia21 def f2_1(tau,faia11,faia21): return ((i_dc + A * math.sin(x))/2 - j - math.sin(faia11/2) - faia21) / BC def f1_2(tau,faia22): return faia22 def f2_2(tau,faia12,faia22): return ((i_dc + A * math.sin(x))/2 - j - math.sin(faia12/2) - faia22) / BC def f1_3(tau,faib2): return faib2 def f2_3(tau,faib1,faib2): return ((i_dc + A * math.sin(x))/2 + j - math.sin(faib1/2)*0.5 - faib2) / BC BC = 0.1 BL = 1 a = 0 b = 2000 N = 200000 h = (b-a)/N s = float() v = float() fai_a = 0 with open("（テスト）φa=0の交流バイアスステップ.txt", "w") as f: for i_dc in np.arange(-1, 1, 0.01): faia11 = 0.0 faia21 = 0.0 faia12 = 0.0 faia22 = 0.0 faib1 = 0.0 faib2 = 0.0 tau = 0.0 for k in range(0, N, 1): tau = k * h x = w * tau j = (faia11 + faia12 - faib1 - 2 * pi * fai_a)/(BL * pi) k1_1 = h * f1_1(tau,faia21) d1_1 = h * f2_1(tau,faia11,faia21) k1_2 = h * f1_2(tau,faia22) d1_2 = h * f2_2(tau,faia12,faia22) k1_3 = h * f1_3(tau,faib2) d1_3 = h * f2_3(tau,faib1,faib2) j = (faia11 + k1_1 + faia12 + k1_2 - (faib1 + k1_3) - 2*pi*fai_a)/(BL*pi) k2_1 = h * f1_1(tau + h/2,faia21 + k1_1/2) d2_1 = h * f2_1(tau + h/2,faia11 + d1_1/2, faia21 + d1_1/2) k2_2 = h * f1_2(tau + h/2,faia22 + k1_2/2) d2_2 = h * f2_2(tau + h/2,faia12 + d1_2/2, faia22 + d1_2/2) k2_3 = h * f1_3(tau + h/2,faib2 + k1_3/2) d2_3 = h * f2_3(tau + h/2,faib1 + d1_3/2, faib2 + d1_3/2) j = (faia11 + k2_1 + faia12 + k2_2 - (faib1 + k2_3) - 2*pi*fai_a)/(BL*pi) k3_1 = h * f1_1(tau + h/2,faia21 + k2_1/2) d3_1 = h * f2_1(tau + h/2,faia11 + d2_1/2, faia21 + d2_1/2) k3_2 = h * f1_2(tau + h/2,faia22 + k2_2/2) d3_2 = h * f2_2(tau + h/2,faia12 + d2_2/2, faia22 + d2_2/2) k3_3 = h * f1_3(tau + h/2,faib2 + k2_3/2) d3_3 = h * f2_3(tau + h/2,faib1 + d2_3/2, faib2 + d2_3/2) j = (faia11 + k3_1 + faia12 + k3_2 - (faib1 + k3_3) - 2*pi*fai_a)/(BL*pi) k4_1 = h * f1_1(tau + h,faia21 + k3_1) d4_1 = h * f2_1(tau + h,faia11 + d3_1, faia21 + d3_1) k4_2 = h * f1_2(tau + h,faia22 + k3_2) d4_2 = h * f2_2(tau + h,faia12 + d3_2, faia22 + d3_2) k4_3 = h * f1_3(tau + h,faib2 + k3_3) d4_3 = h * f2_3(tau + h,faib1 + d3_3, faib2 + d3_3) j = (faia11 + k4_1 + faia12 + k4_2 - (faib1 + k4_3) - 2*pi*fai_a)/(BL*pi) faia11 += 1/6 * (k1_1 + 2 * k2_1 + 2 * k3_1 + k4_1) faia21 += 1/6 * (d1_1 + 2 * d2_1 + 2 * d3_1 + d4_1) faia12 += 1/6 * (k1_2 + 2 * k2_2 + 2 * k3_2 + k4_2) faia22 += 1/6 * (d1_2 + 2 * d2_2 + 2 * d3_2 + d4_2) faib1 += 1/6 * (k1_3 + 2 * k2_3 + 2 * k3_3 + k4_3) faib2 += 1/6 * (d1_3 + 2 * d2_3 + 2 * d3_3 + d4_3) if k > 100000: s += (faia21 + faia22 + faib2)/2 v = s / 100000 print("{:.2f} {:.3f} {:.10f} {:.10f} {:.10f} {:.10f} {:.10f} {:.10f} {:.10f} {:.10f}".format(fai_a, i_dc, v ,j, faia11, faia21, faia12, faia22, faib1, faib2), file=f) s = 0

試したこと

ここに問題に対して試したことを記載してください。

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

現在取り扱っている回路と、正しい結果として得られるはずのグラフ（１）、現在得られている誤ったグラフ（２）を載せます。
グラフの縦軸と横軸が（１）と（２）で入れ替わっているので、ご注意ください。

（（１）の横軸はV、縦軸はIとなっています）
（１）

（２）（２）