Неправильное построение arctg в гнуплоте по точкам

#include <iostream>
#include <cmath>
#include <fstream>
void SLAUSolution(int M, double** A1, double* b1) // 11  A1 = 31x31  B1 = 31
{
    double** A = new double* [M];
    for (int i = 0; i < M; ++i) {
        A[i] = new double[M];
    }
    double* b = new double[M];
    for (int i = 0; i < M; ++i) {
        b[i] = b1[i];
        for (int j = 0; j < M; ++j) {
            A[i][j] = A1[i][j];
        }
    }
    for (int i = 0; i < M - 1; i++) {
        int newi = i;
        for (int k = i + 1; k < M; ++k) {
            if (fabs(A[k][i]) > fabs(A[newi][i])) {
                newi = k;
            }
        }
        for (int k = i; k < M; ++k) {
            double t = A[newi][k];
            A[newi][k] = A[i][k];
            A[i][k] = t;
        }
        {
            double t = b[i];
            b[i] = b[newi];
            b[newi] = t;
        }
        double a0 = A[i][i];
        if (fabs(a0) < 1e-10) {
            std::cout << "Matrix is degenerate" << std::endl;
            return;
        }
        b[i] /= a0;
        for (int j = i; j < M; j++) {
            A[i][j] /= a0;
        }
        for (int k = i + 1; k < M; k++) {
            double ak = A[k][i];
            for (int j = i + 1; j < M; j++) { A[k][j] -= ak * A[i][j]; }
            A[k][i] = 0;
            b[k] -= ak * b[i];
        }
    }
    b[M - 1] /= A[M - 1][M - 1];
    A[M - 1][M - 1] = 1;
    double* xi = new double[M];
    for (int i = M - 1; i >= 0; --i) {
        xi[i] = b[i];
        for (int j = i + 1; j < M; ++j) {
            xi[i] -= A[i][j] * xi[j];
        }
    }
    std::ofstream OUT("slausolution.txt");
    OUT << M << std::endl;
 
    for (int i = 0; i < M; ++i) {
        OUT << xi[i] << std::endl;
    }
    OUT.close();
    // CHECK THE SOLUTION
    double* d;
    d = new double[M];
    for (int i = 0; i < M; ++i) {
        d[i] = -b1[i];
        for (int j = 0; j < M; ++j) {
            d[i] += A1[i][j] * xi[j];
        }
        std::cout << d[i] << std::endl;
    }
    return;
}
void CalculateSplineCoefficients(unsigned N, double* X, double* Y) {
    unsigned N1 = 3 * N - 2; // размер 31
    double** A = new double* [N1];
    double* b = new double[N1]; //размер 31
 
    for (unsigned i = 0; i < N1; i++) // зануляем массивы
    {
        A[i] = new double[N1];
 
        for (unsigned j = 0; j < N1; j++)
        {
            A[i][j] = 0; // 31 на 31
        }
 
        b[i] = 0;
    }
 
    for (unsigned s = 0; s <= N - 2; s++)
    {
        b[s] = Y[s];
        A[s][s] = (X[s + 1] - X[s]) * (X[s + 1] - X[s]) / 6;
        A[s][2 * N - 1 + s] = 1;
    }
 
    for (unsigned s = N - 1; s <= 2 * N - 3; s++)
    {
        unsigned i = s + 2 - N;
        b[s] = Y[i + 1];
        A[s][i + 1] = (X[i + 1] - X[i]) * (X[i + 1] - X[i]) / 6;
        A[s][N + i] = X[i + 1] - X[i];
        A[s][2 * N - 1 + i] = 1;
    }
 
    for (unsigned s = 2 * N - 2; s <= 3 * N - 5; s++)
    {
        unsigned i = s + 3 - 2 * N;
        b[s] = 0;
        A[s][N + i] = 1;
        A[s][N + i - 1] = -1;
        A[s][i] = (X[i - 1] - X[i + 1]) / 2;
    }
 
    A[3 * N - 4][0] = 1;
    A[3 * N - 4][1] = -1;
 
    A[3 * N - 3][N - 2] = 1;
    A[3 * N - 3][N - 1] = -1;
 
    SLAUSolution(N1, A, b);
    std::ifstream INP("slausolution.txt");
    unsigned L;
    INP >> L;
    double* xi = new double[L];
    for (unsigned i = 0; i < L; ++i) {
        INP >> xi[i];
    }
    INP.close();
    std::ofstream OUT("splineresult.txt");
    unsigned K = 10;
    for (unsigned i = 0; i < N - 1; ++i) {
        double Mi = xi[i], Miplus = xi[i + 1], Ci = xi[N + i], di = xi[2 * N - 1 + i];
        for (unsigned k = 0; k < K; ++k) 
        {
            double zk = X[i] * (1 - k * 1 / K) + X[i + 1] * k * 1 / K; //координата х в файле а внизу у
            double fzk = Mi * (X[i + 1] - zk) * (X[i + 1] - zk) * (X[i + 1] - zk) / 6 / (X[i + 1] - X[i]) + Miplus * (zk - X[i]) * (zk - X[i]) * (zk - X[i]) / 6 / (X[i + 1] - X[i]) + Ci * (zk - X[i]) + di;
            OUT << zk << " " << fzk << std::endl;
        }
    }
    OUT.close();
    return;
}
 
void Test(const char* fname) {
    std::ofstream OFOF(fname);
    unsigned N = 11;
    OFOF << N << std::endl;
    for (double i = 0; i < N; i += 1) {
        OFOF << i << " " << std::atan(i) << std::endl;
    }
    OFOF.close();
    return;
}
 
int main()
{
    unsigned N;
    const char fname[] = "Krik";
    Test(fname);
    std::ifstream AFAF(fname);
    AFAF >> N;
    double* X = new double[N];
    double* Y = new double[N];
    for (unsigned i = 0; i < N; ++i) 
    {
        AFAF >> X[i] >> Y[i];
    }
    CalculateSplineCoefficients(N, X, Y);
    return 0;
    double** A = new double* [2];
    A[0] = new double[2];
    A[1] = new double[2];
    double* b = new double[2];
    A[0][0] = 0;
    A[0][1] = 1;
    A[1][0] = 1;
    A[1][1] = 0;
    b[0] = 1;
    b[1] = -1;
    SLAUSolution(2, A, b);
    return 0;
}