# Trying to implement Newton's interpolation polynomial - not evaluating divided differences

I'm trying to find Newton's interpolation polynomial for a given function and a set of points. However, I'm having some issues with finding the divided differences (in this case I've notated them with R) with nested Do expressions and iterating over the difference in indices between the first and last point. Whenever I run it for x != y, simply R[x, y] is printed.

n = 6;
(* Given function *)
f[x_] := 1/(1 + 4 x);
(* Given points *)
Do[x[k] = 0.5*k, {k, 0, n}];

(* f[x_i] = f(x_i) *)
Do[R[i_, i_] := f[x[i]], {i, 0, n}];

(* f[x_0, ... , x_n] = *)
(* (f[x_1, ... , x_n] - f[x_0, ... , x_n-1]) / (x_n - x_0) *)
Do[
Do[R[i_, i_ + d_] = (R[i + 1, i + d] - R[i, i + d - 1]) / (x[i + d] - x[i]),
{i, 0, n - d}],
{d, 1, n}];

Prod[t_, 0] = 1;
Do[Prod[t_, k_] := Product[(t - x[l]), {l, 0, k - 1}], {k, 0, n}];

L[f_, t_] := Sum[R[0, k]*Prod[t, k], {k, 0, n}];

(R[1, 1] - R[0, 0])/(x[1] - x[0])
R[0, 1]


And yes, I know there is a function for interpolating polynomials, I'm trying to accomplish finding the polynomial without it.

This is not a direct answer for your problem, nor the most efficient approach, but...

This shows how easily the idea of Newton divided differences may be expressed in Mathematica (the formula was lifted from your code-comment, which inspired this approach).

f[x_] := 1/(1 + 4 x);
InternalInheritedBlock[{f},
f[x0_, xx___, xn_] := (f[xx, xn] - f[x0, xx])/(xn - x0);
f @@ (0.5 Range[0, 6])
]
(*  0.03031043031043034 *)


Check with InterpolatingPolynomial:

InterpolatingPolynomial[Transpose[{#, f[#]}] &[0.5 Range[0, 6]], x];
Level[%, {2*6}]
(*  {0.03031043031043036, -1. + x}  *)