Newton’s Forward Interpolation

Write a program in C to find 𝒚 at 𝒙 = 𝟏. 𝟓 using the Newton’s Forward Difference Method given data:
x: 1, 2, 3, 4
y: 1, 8, 27, 64

#include <stdio.h>
#include<conio.h>
int main() {
 clrscr();
 double y[4] = {1, 8, 27, 64}, x0 = 1, h = 1, xp = 1.5, u, yp;
 double d[4][4];
 int i, j;
 for (i = 0; i < 4; i++) d[i][0] = y[i];
 for (j = 1; j < 4; j++)
 for (i = 0; i < 4 - j; i++)
 d[i][j] = d[i + 1][j - 1] - d[i][j - 1];
 u = (xp - x0) / h;
 yp = d[0][0]
+ u * d[0][1]
+ u * (u - 1) * d[0][2] / 2.0
+ u * (u - 1) * (u - 2) * d[0][3] / 6.0;
 printf("\n\n\ty(%.1f) = %.3f\n", xp, yp);
 getch();
 return 0;
}

Output