I’m trying to fix my Hermite interpolation program, which is supposed to output the coefficients and value at a point, but it seems the results are wrong. It currently looks like this:

#include <iostream>
#include <vector>

using namespace std;

struct DataPoint {
    double x;
    vector<double> values;  
};

void computeDividedDifferences(vector<DataPoint>& points, vector<int>& multiplicity, vector<vector<double>>& divDiff) {
    int n = points.size();
    int totalNodes = 0;
    for (int i = 0; i < n; ++i) {
        totalNodes += multiplicity[i];
    }

    for (int i = 0; i < totalNodes; ++i) {
        for (int j = 0; j < totalNodes; ++j) {
            divDiff[i][j] = 0.0;
        }
    }

    int index = 0;
    for (int i = 0; i < n; ++i) {
        for (int k = 0; k < multiplicity[i]; ++k) {
            divDiff[index][0] = points[i].values[k];
            if (k != 0) {
                divDiff[index][1] = points[i].values[k] / k;
            } else if (index != 0) {
                divDiff[index][1] = (divDiff[index][0] - divDiff[index - 1][0]) / (points[i].x - points[(index - 1) / multiplicity[(index - 1) / multiplicity[i]]].x);
            }
            index++;
        }
    }

    for (int j = 2; j < totalNodes; ++j) {
        for (int i = j; i < totalNodes; ++i) {
            divDiff[i][j] = (divDiff[i][j - 1] - divDiff[i - 1][j - 1]) / (points[i / multiplicity[i / multiplicity[i]]].x - points[(i - j) / multiplicity[(i - j) / multiplicity[i]]].x);
        }
    }
}

double hermitePolynomial(double x, vector<DataPoint>& points, vector<int>& multiplicity, vector<vector<double>>& divDiff) {
    int n = points.size();
    int totalNodes = 0;
    for (int i = 0; i < n; ++i) {
        totalNodes += multiplicity[i];
    }

    double result = divDiff[0][0];
    double product = 1.0;

    for (int i = 1; i < totalNodes; ++i) {
        product *= (x - points[i / multiplicity[i / multiplicity[i]]].x);
        result += divDiff[i][i] * product;
    }

    return result;
}

void printCoefficients(vector<DataPoint>& points, vector<int>& multiplicity, vector<vector<double>>& divDiff) {
    int n = points.size();
    int totalNodes = 0;
    for (int i = 0; i < n; ++i) {
        totalNodes += multiplicity[i];
    }

    cout << "Hermite polynomial coefficients:" << endl;
    for (int i = 0; i < totalNodes; ++i) {
        cout << "a" << i << " = " << divDiff[i][i] << endl;
    }
}

int main() {
    int k;
    cout << "Enter the number of nodes minus one (k): ";
    cin >> k;
    int n = k + 1;

    vector<DataPoint> points(n);
    vector<int> multiplicity(n);

    cout << "Enter the values of interpolation nodes (x[i], where i = 0, ..., k):" << endl;
    for (int i = 0; i < n; ++i) {
        cout << "x[" << i << "] = ";
        cin >> points[i].x;
    }

    cout << "Enter the multiplicities of nodes (m[i], where i = 0, ..., k):" << endl;
    for (int i = 0; i < n; ++i) {
        cout << "m[" << i << "] = ";
        cin >> multiplicity[i];
    }

    for (int i = 0; i < n; ++i) {
        points[i].values.resize(multiplicity[i]);
        cout << "Enter the function values and its derivatives at node " << points[i].x << " (f[i][s], where s = 0, ..., m[i] - 1):" << endl;
        for (int j = 0; j < multiplicity[i]; ++j) {
            cout << "f[" << i << "][" << j << "] = ";
            cin >> points[i].values[j];
        }
    }

    vector<vector<double>> divDiff(2 * n, vector<double>(2 * n));
    computeDividedDifferences(points, multiplicity, divDiff);

    double x;
    cout << "Enter the value of x to calculate the Hermite polynomial: ";
    cin >> x;

    double result = hermitePolynomial(x, points, multiplicity, divDiff);
    cout << "The value of the Hermite polynomial at point " << x << " is: " << result << endl;

    printCoefficients(points, multiplicity, divDiff);

    return 0;
}

The input i give is:

Enter the number of nodes minus one (k): 2
Enter the values of interpolation nodes (x[i], where i = 0, ..., k):
x[0] = 0
x[1] = 1
x[2] = 2
Enter the multiplicities of nodes (m[i], where i = 0, ..., k):
m[0] = 1
m[1] = 1
m[2] = 2
Enter the function values and its derivatives at node 0 (f[i][s], where s = 0, ..., m[i] - 1):
f[0][0] = 2
Enter the function values and its derivatives at node 1 (f[i][s], where s = 0, ..., m[i] - 1):
f[1][0] = 3
Enter the function values and its derivatives at node 2 (f[i][s], where s = 0, ..., m[i] - 1):
f[2][0] = 4
f[2][1] = 5
Enter the value of x to calculate the Hermite polynomial: 4

the results i get are:

The value of the Hermite polynomial at point 4 is: 3.79421e+99
Hermite polynomial coefficients:
a0 = 2
a1 = 1
a2 = 0
a3 = 1.58092e+98

But the results that another Hermite interpolation calculator I’m trying to recreate gives are

a[0] =  2 
a[1] = -5 
a[2] =  6 
a[3] = -2 
The value of polynomial at  x =  4: 
H(x) =  5.4

I tried to debug it, read more about how Hermite interpolation works and look at the code again, but I’m still not sure what seems to be the problem. I want the code to output correct coefficients and value at a point.