function [x,y] = myFD(alpha, beta, a, b, N) % 求解常微分方程边值问题 y'' + p(x)y' + q(x)y = r(x), a < x < b % y(a) = alpha, y(b) = beta % N 表示对区间 [a,b] 剖分 N 等份 % 返回值 x 为 [a,b] 的剖分节点 % 返回值 y 为剖分节点上的数值解 %p(x)=3 q(x)=2 y(x)=sin(x)

首先，我们需要将二阶常微分方程转化为一阶常微分方程组。设 $y_1(x)=y(x)$，$y_2(x)=y'(x)$，则方程组可以表示为：

$$ \begin{cases} y_1'(x) = y_2(x) \\ y_2'(x) = r(x) - p(x)y_2(x) - q(x)y_1(x) \end{cases} $$

其中，$y_1(a)=\alpha$，$y_1(b)=\beta$。

接下来，我们采用有限差分方法（Finite Difference Method）来数值求解上述方程组。将区间 $[a,b]$ 等分成 $N+1$ 个节点，其中 $a=x_0<x_1<\cdots<x_N=b$，$h=\frac{b-a}{N}$ 为步长。对于节点 $x_i$，我们用 $y_{1,i}$ 和 $y_{2,i}$ 分别表示 $y_1(x_i)$ 和 $y_2(x_i)$ 的近似值。

根据二阶导数的定义，我们可以使用中心差分公式来近似求解 $y_2'(x_i)$：

$$y_2'(x_i)\approx\frac{y_{2,i+1}-y_{2,i-1}}{2h}$$

同样地，我们可以使用中心差分公式来近似求解 $y_1'(x_i)$：

$$y_1'(x_i)\approx\frac{y_{1,i+1}-y_{1,i-1}}{2h}$$

将上述两个式子带入到方程组中，我们可以得到：

$$ \begin{cases} y_{1,i+1}=y_{1,i}+hy_{2,i}+\frac{h^2}{2}\left(r(x_i)-p(x_i)y_{2,i}-q(x_i)y_{1,i}\right) \\ y_{2,i+1}=y_{2,i}+\frac{h}{2}\left(r(x_i)-p(x_i)y_{2,i}-q(x_i)y_{1,i}\right)+\frac{h}{2}\left(r(x_{i+1})-p(x_{i+1})y_{2,i+1}-q(x_{i+1})y_{1,i+1}\right) \end{cases} $$

注意到，在 $y_{2,i+1}$ 的计算中，我们需要使用 $y_{1,i+1}$ 的值。因此，我们需要使用一个迭代方法来不断地更新 $y_{1,i+1}$ 和 $y_{2,i+1}$ 直到收敛。这里我们采用简单的 Jacobi 迭代法，即每次更新时，使用上一次迭代得到的 $y_{1,i+1}^{(k)}$ 和 $y_{2,i+1}^{(k)}$ 来计算出新的 $y_{1,i+1}^{(k+1)}$ 和 $y_{2,i+1}^{(k+1)}$。具体来说，迭代公式如下：

$$ \begin{cases} y_{1,i+1}^{(k+1)}=y_{1,i}+hy_{2,i}+\frac{h^2}{2}\left(r(x_i)-p(x_i)y_{2,i}-q(x_i)y_{1,i}^{(k)}\right) \\ y_{2,i+1}^{(k+1)}=y_{2,i}+\frac{h}{2}\left(r(x_i)-p(x_i)y_{2,i}-q(x_i)y_{1,i}^{(k)}\right)+\frac{h}{2}\left(r(x_{i+1})-p(x_{i+1})y_{2,i+1}^{(k)}-q(x_{i+1})y_{1,i+1}^{(k+1)}\right) \end{cases} $$

其中，$k$ 表示迭代次数。

最后，我们可以通过以上迭代方法，依次计算出 $y_{1,1},y_{2,1},y_{1,2},y_{2,2},\cdots,y_{1,N},y_{2,N}$ 的近似值，从而得到整个区间 $[a,b]$ 上的数值解 $y(x)$。代码如下：