多次元の非線形連立方程式の解き方

Question

前提・実現したいこと

現在、非線形波動方程式の数値計算をしています。
その際に多次元の非線形連立方程式を解かなければいけません。
そこでpython3のfsolveを使ってコードを書きました。

取り組みの手順として、
まず多次元の線形連立方程式をfsolveで解き、非線形連立方程式を解こうと考えました。
それぞれのコードを書き実行した結果次のようなエラー・警告・問題点が出ました。

・[線形連立方程式]
１．RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last ten iterations.
warnings.warn(msg, RuntimeWarning)
２．出力値が違う。

・[非線形連立方程式]
TypeError: 'int' object is not callable

・[できれば、、、]
できれば、より早く計算を行いたいと考えています。
※方程式を解けていないため検討をしていません。
しかし、疎行列を扱うためscipy.sparseを用いればよいのではないかと考えています。

それぞれを解決したいです。
下の非線形連立方程式を解ければよいのでfsolve以外でも大丈夫です。
どうぞよろしくお願い致します。
###検討する非線形連立方程式
手書きで見にくいですが解きたいスキームとそれを行列で表した式です。
m:時間空間を分割した際のステップ数
k:x空間を分割した値

該当のソースコード

%matplotlib nbagg
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy import optimize
from scipy.linalg import toeplitz
from scipy.optimize import fsolve

T=20
N=2000
dt=T/N
L=6
K=300
dx=L/K
r=dt**2/(dx**2)
R=dt**2/dx
u=np.zeros((K+1,N,))
#二階差分
D2=toeplitz([-2, 1, *np.zeros(K-1)])

A=2/r*np.eye(K+1)-D2
A[0,0]=1+1/R
A[0,1]=1/R-1
A[K,K]=1/R+1
A[K,K-1]=1/R-1
B=4/r*np.eye(K+1)
B[0,0]=-1/R
B[0,1]=-1/R
B[K,K]=-1/R
B[K,K-1]=-1/R
C=-2/r*np.eye(K+1)+D2
C[0,0]=-1-1/R
C[0,1]=1-1/R
C[K,K]=-1-1/R
C[K,K-1]=1-1/R

#初期値の設定
for k in range(0,K+1):
    u[k,0]=np.exp(-(k*dx-L/2)**2)
    u[k,1]=u[k,0]-dt*2*(k*dx-L/2)**2*np.exp(-(k*dx-L/2)**2)

#----------------------非線形（線形）連立方程式を解く-------------------
def func(x):
    x1=u[:,n-1]
    x2=u[:,n]
    x3=u[:,n+1]
#非線形の場合-1/4(x1**3+x3**2*x1+x3*x1**2+x1**3)をつける
    return  np.dot(A,x3)-np.dot(B,x2)-np.dot(C,x1)#-1/4(x1**3+x3**2*x1+x3*x1**2+x1**3)
#数値計算
for n in range(1,N-1):
    x0=u[:,n]
    u[:,n+1]=fsolve(func, x0)
#--------------------------------------------------------------------

# matplotlibで表示
x = np.linspace(0, L, K+1)  # maplotlibで表示するためのx軸の設定
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(1,1,1)
# アニメ更新用の関数
def update_func(i):
    # 前のフレームで描画されたグラフを消去
    ax.clear()
    ax.plot(x, u[:,i], color='red')#力学的境界条件
    # 軸の設定
    ax.set_ylim(-1.2, 1.2)
    # 軸ラベルの設定
    ax.set_xlabel('x', fontsize=12)
    ax.set_ylabel('u', fontsize=12)
ani = animation.FuncAnimation(fig, update_func, frames=N, interval=1,repeat=False)
#ani.save("output_riki_3.mp4", writer="ffmpeg", fps=30, bitrate=1000)
plt.show()

試したこと

１．別の手法を使って線形連立方程式を解き、解があっているかどうかを確かめました。
そのため、上のコードの「非線形（線形）連立方程式をとく」の中身がエラー・警告の原因だと考えています。
下記に線形連立方程式を解いたコードを書きます。
２．インターネット上では3次元の非線形連立方程式をとくサイトはありましたが多次元のサイトはなく、見よう見まねで先ほど書いたコードを書きました。

%matplotlib nbagg
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy import optimize
from scipy.linalg import toeplitz
import scipy.sparse
import scipy.sparse.linalg as spla

T=20
N=2000
dt=T/N
L=6
K=300
dx=L/K
r=dt**2/(dx**2)
R=dt**2/dx
u=np.zeros((K+1,N,))

A=2/r*np.eye(K+1)-D2
A[0,0]=1+1/R
A[0,1]=1/R-1
A[K,K]=1/R+1
A[K,K-1]=1/R-1
B=4/r*np.eye(K+1)
B[0,0]=-1/R
B[0,1]=-1/R
B[K,K]=-1/R
B[K,K-1]=-1/R
C=-2/r*np.eye(K+1)+D2
C[0,0]=-1-1/R
C[0,1]=1-1/R
C[K,K]=-1-1/R
C[K,K-1]=1-1/R

#初期値の設定
for k in range(0,K+1):
    u[k,0]=np.exp(-(k*dx-L/2)**2)
    u[k,1]=u[k,0]-dt*2*(k*dx-L/2)*np.exp(-(k*dx-L/2)**2)

#数値計算
for n in range(1,N-1):
    b=np.dot(B,u[:,n])+np.dot(C,u[:,n-1])
    a_csr = scipy.sparse.csr_matrix(A)
    u[:,n+1] = spla.dsolve.spsolve(a_csr, b)

# matplotlibで表示
x = np.linspace(0, L, K+1)  # maplotlibで表示するためのx軸の設定
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(1,1,1)
# アニメ更新用の関数
def update_func(i):
    # 前のフレームで描画されたグラフを消去
    ax.clear()
    ax.plot(x, u[:,i], color='red')
    #軸の設定
    ax.set_ylim(-1.2, 1.2)
    # 軸ラベルの設定
    ax.set_xlabel('x', fontsize=12)
    ax.set_ylabel('u', fontsize=12)
ani = animation.FuncAnimation(fig, update_func, frames=N, interval=1,repeat=False)
#ani.save("output_riki_3.mp4", writer="ffmpeg", fps=30, bitrate=1000)
plt.show()