Реализация интерполяционной формулы Ньютона

#include <iostream>
#include <conio.h>
 
using namespace std;
 
int main() {
    const int ixy = 21;
    int ixy0 = 21;
 
    double x[ixy], ch = 1.0;
 
    double xs[2];
    cout << "Enter x of search 1: ";
    cin >> xs[0];
    cout << "Enter x of search 2: ";
    cin >> xs[1];
 
    for (int i = 0; ch <= 3.0, i < ixy0; i++, ch += 0.1)
        x[i] = ch;
 
    double y[] = {
        1.2661, 1.3262, 1.3937, 1.4693, 1.5535,
        1.6467, 1.7500, 1.8640, 1.9896, 2.1277,
        2.2796, 2.4463, 2.6291, 2.8296, 3.0493,
        3.2898, 3.5533, 3.8417, 4.1573, 4.5027,
        4.8808
    };
 
    double  d[100][100];
    for (int i = 0; i < 21; i++){
        d[i][0] = y[i];
    }
 
    int coutm = 21;
    int coutp = 0;
 
    for (int p = 1; p < 5; p++){
        coutm -= 1;
        coutp += 1;
        for (int i = 0; i < coutm; i++){
            d[i][p] = d[i + 1][p - 1] - d[i][p - 1];
        }
    }
    
    int cs = 0;
 
    for (int j = 0; j < 21; j++) {
        if (x[j] == xs[0])
            break;
        cout << "x[j]: " << x[j] << "\t" << xs[0] << "\n" << cs << "\n\n";
        cs++;
    }
 
 
    double rez = 0;
 
    for (int k = 1; k <= coutp; k++){
        rez += pow(-1, k + 1) * d[cs][k] / k;
    }
 
    cout << rez;
 
    int cst = 0;
 
    _getch();
 
    return 0;
}