Интерполяция функции кубическим сплайном

#include <iostream>
#include <iomanip>
 
using namespace std;
 
void CubicSpline(double a, double b, double dx, int n)
{
    //Вычисление шага
    double h = (b - a) / n;
 
    //Число xj на интервале [a,b]
    int count = (int)((b - a) / dx + 1);
    cout << endl << h << endl << count << endl << endl;
 
    double* x = new double[n + 1];
    double* y = new double[n + 1];
    cout << endl;
    cout << "x[i]" << "\t" << "F[x[i]]" << endl;
    for (int i = 0; i <= n; i++)
    {
        x[i] = a + (double)i * h;
        if (x[i] == 0) y[i] = 0;
        else y[i] = sin(pow(x[i], 4)) / (pow(x[i], 2));
        cout << fixed << setprecision(17) << x[i] << "\t" << y[i] << endl;
    }
 
    cout << endl << "NAKLONI" << endl;
    //Массив значений наклонов в точках
    double* m = new double[n + 1];
    m[0] = (4.0 * y[1] - y[2] - 3.0 * y[0]) / 2.0 * h;
    m[n] = (3.0 * y[n] + y[n - 2] - 4.0 * y[n - 1]) / (2.0 * h);
    for (int i = 1; i < n; i++)
    {
        m[i] = (y[i + 1] - y[i - 1]) / (2.0 * h);
    }
    cout << m[0] << endl;
    for (int i = 1; i < n; i++)
    {
        cout << fixed << setprecision(17) << m[i] << endl;
    }
 
    cout << endl;
 
 
    //Нахождение точек интерполянта кубическим сплайном
    double* spline = new double[(int)count];
    double firstTerm(0), secondTerm(0), thirdTerm(0), fourthTerm(0);
    int i(0), j(0);
 
    FILE* OutFile;
    fopen_s(&OutFile, "res_interp.txt", "wt");
 
    // чтобы очистить файл
    FILE *file_clear;
    fopen_s(&file_clear, "res.txt", "wt");
    fprintf(file_clear, "");
    fclose(file_clear);
 
    printf("\n\n\n%s\t\t%s\t\t%s\n", "xj", "F(xj)", "f(xj)-F(xj)");
    for (int i = 0; i <= count; i++)
    {
        double xj = a + i * dx;
 
 
        if (xj >= x[i] && xj <= x[i + 1])
 
 
        {
 
            //Вычисление кубического сплайна
            firstTerm = (pow((x[i + 1] - xj), 2) * (2 * (xj - x[i]) + h) * y[i]) / (h * h * h);
            secondTerm = (pow((xj - x[i]), 2) * (2 * (x[i + 1] - xj) + h) * y[i + 1]) / (h * h * h);
            thirdTerm = (pow((x[i + 1] - xj), 2) * (xj - x[i]) * m[i]) / (h * h);
            fourthTerm = (pow((xj - x[i]), 2) * (xj - x[i + 1]) * m[i + 1]) / (h * h);
 
            spline[j] = firstTerm + secondTerm + thirdTerm + fourthTerm;
 
 
 
 
 
 
 
            FILE *file_out;
            fopen_s(&file_out, "res.txt", "a+");
            fclose(file_out);
        }
    }
    fclose(OutFile);
 
 
 
 
 
    for (int i = 0; i <= count; i++)
    {
        double xj = a + i * dx;
        // подбор k, где s в интервале [x_k-1; x_k]
 
 
        if (xj >= x[i] && xj <= x[i + 1])
 
 
        {
 
            //Вычисление кубического сплайна
            firstTerm = (pow((x[i + 1] - xj), 2) * (2 * (xj - x[i]) + h) * y[i]) / (h * h * h);
            secondTerm = (pow((xj - x[i]), 2) * (2 * (x[i + 1] - xj) + h) * y[i + 1]) / (h * h * h);
            thirdTerm = (pow((x[i + 1] - xj), 2) * (xj - x[i]) * m[i]) / (h * h);
            fourthTerm = (pow((xj - x[i]), 2) * (xj - x[i + 1]) * m[i + 1]) / (h * h);
 
            spline[j] = firstTerm + secondTerm + thirdTerm + fourthTerm;
 
            cout << xj << "\t" << spline[j] << "\t" << y[i] - spline[j] << endl;
        }
 
    }
 
 
 
 
    delete[] x;
    delete[] y;
    delete[] m;
    delete[] spline;
 
}
 
 
void main()
{
 
    //Интервал
    double a(0), b(1.68);
    //Число узлов
    int n(11);
    //Xj c шагом dx
    double dx = 0.08;
 
    CubicSpline(a, b, dx, n);
    system("pause");
}