差分法求解微分方程

import numpy as np
from numpy import exp
from sympy import symbols, solve, Eq
from math import e
from matplotlib import pyplot as plt

# d2phi/dx2 + d phi/dx - 2* phi =0, phi(0) = 0, phi(4)=1

# 通用部分开始
x_l, p_x_l = 0, 0  # 左端 边界条件
x_r, p_x_r = 4, 1  # 右端 边界条件

n = 17  # 节点数, n>=3
X = np.linspace(x_l, x_r, n)
dx = X[1] - X[0]
P = [symbols(f'p{i}') for i in range(n)]  # 定义 p0, p1, ..., pn # p for phi
Dp_dx = [symbols(f'dp_dx{i}') for i in range(n)]  # 一阶导数
D2p_dx2 = [symbols(f'd2p_dx2{i}') for i in range(n)]  # 二阶导数
# print(P)
# print(X)

#  左边界处
Dp_dx[0] = (-3 * P[0] + 4 * P[1] - P[2]) / (2 * dx)
D2p_dx2[0] = (P[0] - 2 * P[1] + P[2]) / (dx * dx)

# 中间节点
for i in range(1, n-1):
    Dp_dx[i] = (P[i+1] - P[i-1]) / (2 * dx)
    D2p_dx2[i] = (P[i+1] - 2*P[i] + P[i-1]) / (dx * dx)

# 右边界
Dp_dx[n-1] = (3*P[n-1] - 4*P[n-2] + P[n-3]) / (2 * dx)
D2p_dx2[n-1] = (P[n-1] - 2 * P[n-2] + P[n-3]) / (dx * dx)
# 通用部分结束
eqs = [Eq(P[0], p_x_l), Eq(P[n-1], p_x_r)]  # 边界条件
for i in range(1, n-1):  # 内部满足微分方程
    eq = Eq(D2p_dx2[i] + Dp_dx[i] - 2 * P[i], 0)
    eqs.append(eq)
# print(eqs)
r = solve(eqs,P)
# print(r)

p = np.zeros(n)
for i in range(n):
    p[i] = r[P[i]]
print(p)


x_ = np.linspace(x_l, x_r,1000)
y = (np.exp(x_) - np.exp(-2*x_)) / (e**4 - e**(-8))  # 真实解
plt.plot(x_, y, c="g", label="真实解")

x = np.linspace(x_l, x_r, n)
plt.scatter(x, p, c="r", label="差分解")
plt.legend(loc="lower right")
plt.show()