对于一个如下的一元二阶常微分方程
$$\{\begin{array}{l}-u^{\prime \prime }+q\left(x\right)u=f\left(x\right)a<x<b\\ u\left(a\right)=\alpha u\left(b\right)=\beta \end{array}$$
其中二阶导数可以换成差分格式
$$u^{\prime \prime }\left(x_i\right)=\frac{1}{h^2}\left(u\left(x_{i-1}\right)-2u\left(x_i\right)+u\left(x_{i+1}\right)\right)$$
将$u(x_i)$简写为$u_i$
$$-\frac{u_{i-1}}{h^2}+\left(q\left(x_i\right)+\frac{2}{h^2}\right)u_i-\frac{u_{i+1}}{h^2}=f\left(x_i\right)$$
将a到b的区域平均分成m份，解以上的式子可以看做解一个三对角方程组
$$\left(\begin{array}{cc}\begin{array}{ccc}2+h^{2}q(x_1)& -1& \\ -1& 2+h^{2}q(x_2)& -1\\ & \ddots & \ddots \end{array}& \begin{array}{ccc}& & \\ & & \\ \ddots & & \end{array}\\ \begin{array}{ccc}& & \\ & & \\ & & \end{array}& \begin{array}{ccc}-1& 2+h^{2}q(x_{m-2})& -1\\ & -1& 2+h^{2}q(x_{m-1})\\ & & \end{array}\end{array}\right)\left(\begin{array}{c}{u}_1\\ u_2\\ \begin{array}{c}⋮\\ u_{m-2}\end{array}\\ u_{m-1}\end{array}\right)=\left(\begin{array}{c}{h}^{2}f\left(x_1\right)+\alpha \\ h^{2}f\left(x_2\right)\\ \begin{array}{c}⋮\\ h^{2}f\left(x_{m-2}\right)\end{array}\\ h^{2}f\left(x_{m-1}\right)+\beta \end{array}\right)$$
这里是左右两边同时乘以$h^2$之后展开成的矩阵形式，取$f(x)=e^x(\sin (x)-2\cos (x))$ $q(x)=1$代码如下

#include<stdio.h>
#include<math.h>
#define M_PI 3.14159265358979323846

#define M 20 //m等分
#define N M-1
/*
a、b、c是三对角矩阵的下对角线、主对角线和上对角线的元素，d是右侧向量，x是结果向量
*/
void thomas(double a[], double b[], double c[], double d[], double x[])
{
    double y[N],p[N],q[N];
    int i;
    p[0]=b[0];
    for(i=0;i<N-1;i++){
        q[i]=c[i]/p[i];
        p[i+1]=b[i+1]-a[i]*q[i];
    }
    y[0]=d[0]/p[0];
    for(i=1;i<N;i++){
        y[i]=(d[i]-a[i-1]*y[i-1])/p[i];
    }
    x[N-1]=y[N-1];
    for(i=N-2;i>=0;i--){
        x[i]=y[i]-q[i]*x[i+1];
    }
}

double f(double x){
    return exp(x)*(sin(x)-2*cos(x));
}

double q(double x){
    return 1;
}
int main(){
    double a=0.0;//下限
    double b=M_PI;//上限
    double h=(b-a)/M;//步长
    double alpha=0.0;
    double beta=0.0;
    double zhu[N],xia[N-1],shang[N-1],you[N],x[N];
    int i;
    for(i=0;i<N;i++){
        zhu[i]=q(a+i*h+h)*h*h+2;
        you[i]=f(a+i*h+h)*h*h;
    }
    you[0]+=alpha;
    you[N-1]+=beta;
    for(i=0;i<N-1;i++){
        xia[i]=-1;
        shang[i]=-1;
    }
    thomas(xia,zhu,shang,you,x);
    for(i=0;i<N;i++){
        printf("u[%d]=%lf\n",i+1,x[i]);
    }
    return 0;
}

这里也包含了追赶法求三对角矩阵的函数