请用C++实现matlab的spline函数并且编写一个例程

以下是一个使用C++实现matlab的spline函数的例程。

#include <iostream>
#include <vector>
#include <map>

using namespace std;

struct SplineCoef {
    double a, b, c, d;
};

vector<SplineCoef> spline(vector<double>& x, vector<double>& y) {
    int n = x.size();
    vector<double> h(n - 1), alpha(n - 1), l(n), mu(n - 1), z(n);
    vector<SplineCoef> coef(n - 1);

    for (int i = 0; i < n - 1; i++) {
        h[i] = x[i + 1] - x[i];
        alpha[i] = 3.0 / h[i] * (y[i + 1] - y[i]) - 3.0 / h[i - 1] * (y[i] - y[i - 1]);
    }

    l[0] = 1.0;
    mu[0] = 0.0;
    z[0] = 0.0;

    for (int i = 1; i < n - 1; i++) {
        l[i] = 2.0 * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1];
        mu[i] = h[i] / l[i];
        z[i] = (alpha[i - 1] - h[i - 1] * z[i - 1]) / l[i];
    }

    l[n - 1] = 1.0;
    z[n - 1] = 0.0;
    coef[n - 2].c = 0.0;

    for (int j = n - 3; j >= 0; j--) {
        coef[j].c = z[j] - mu[j] * coef[j + 1].c;
        coef[j].b = (y[j + 1] - y[j]) / h[j] - h[j] * (coef[j + 1].c + 2.0 * coef[j].c) / 3.0;
        coef[j].d = (coef[j + 1].c - coef[j].c) / (3.0 * h[j]);
        coef[j].a = y[j];
    }

    return coef;
}

double spline_eval(vector<SplineCoef>& coef, double x) {
    int n = coef.size();
    if (x <= coef[0].a) return coef[0].a;
    if (x >= coef[n - 1].a) return coef[n - 1].a;

    int i = 0, j = n - 1;
    while (i < j) {
        int k = (i + j) / 2;
        if (x <= coef[k].a) j = k;
        else i = k + 1;
    }

    double dx = x - coef[i].a;
    return coef[i].a + coef[i].b * dx + coef[i].c * dx * dx + coef[i].d * dx * dx * dx;
}

int main() {
    vector<double> x = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0};
    vector<double> y = {0.0, 0.8, 0.9, 0.1, -0.8, -1.0};

    vector<SplineCoef> coef = spline(x, y);

    // 输出样条插值系数
    for (int i = 0; i < coef.size(); i++) {
        cout << "a = " << coef[i].a << ", b = " << coef[i].b << ", c = " << coef[i].c << ", d = " << coef[i].d << endl;
    }

    // 输出样条插值结果
    for (double i = 0.0; i <= 5.0; i += 0.1) {
        cout << "x = " << i << ", y = " << spline_eval(coef, i) << endl;
    }

    return 0;
}

该例程先定义了样本点的x和y坐标，然后调用spline函数得到样条插值系数，最后使用spline_eval函数得到样条插值结果。