Off[General::spell1];
Clear[DividedDifference];
DividedDifferece[XY0_] :=
Block[{j,k,n,a,p,sum,X,Y},
  X = Transpose[XY0][[1]];
  Y = Transpose[XY0][[2]];
  n = Length[XY0]-1;
  d = Table[0,{n+1},{n+1}];
  a = Table[0,{n+1}];
  p = Table[0,{n+1}];
  For[k=0,k<=n,k++,
    d[[k+1,1]] = Y[[k+1]];];
  For[j=1,j<=n,j++,
    For[k=j,k<=n,k++,
      d[[k+1,j+1]]=
       (d[[k+1,j]]- d[[k,j]])/
         (X[[k+1]]-X[[k-j+1]]) ];];];
  On[General::spell1];

Clear[data];
data={{0,-1}, {1, 2}, {2,1}, {3, 0}}

N[DividedDifferece[data]];
TableForm[d]

Clear[data];
data={{0,1}, {0.5, 2}, {1,0}, {1.5, -1}, {2, 1}}

N[DividedDifferece[data]];
TableForm[d]

Clear[f,x];
f[x_] = Cos[x]
xy = Table[{i,N[f[i],8]},{i,0,3}]

N[DividedDifferece[xy]]; 
TableForm[d]

Off[General::spell1];
Clear[NewtonForm];
NewtonForm[XY0_] :=
Block[{j,k,n,a,p,sum,X,Y},
  X = Transpose[XY0][[1]];
  Y = Transpose[XY0][[2]];
  n = Length[XY0]-1;
  d = Table[0,{n+1},{n+1}];
  a = Table[0,{n+1}];
  p = Table[0,{n+1}];
  For[k=0,k<=n,k++,
    d[[k+1,1]] = Y[[k+1]];];
  For[j=1,j<=n,j++,
    For[k=j,k<=n,k++,
      d[[k+1,j+1]]=
      (d[[k+1,j]]-d[[k,j]])/
           (X[[k+1]]-X[[k-j+1]]) ];];
  For[k=0,k<=n,k++,
    a[[k+1]] = Together[d[[k+1,k+1]]];];
  p[[1]] = 1;
  For[k=1,k<=n,k++,
    p[[k+1]] = p[[k]] (x-X[[k]]) ];
  sum = 0;
  For[k=0,k<=n,k++,
    sum = sum + a[[k+1]] p[[k+1]];];
  sum = sum/.h->H;
  Return[sum]; ];
On[General::spell1];

Clear[data];
data={{0,-1}, {1, 2}, {2,1}, {3, 0}}

N[DividedDifferece[data]];
TableForm[d]

Clear[P,x];
P[x_] = N[NewtonForm[data],10]

dots = ListPlot[data,PlotStyle->{PointSize[0.02]},
        DisplayFunction->Identity];
g = Plot[{P[x]}, {x, -0.2, 3.5},
DisplayFunction->Identity];
Show[g,dots, PlotLabel->"data and the Newton polynomial.",
AxesLabel->{"x","y"},DisplayFunction->$DisplayFunction];

Clear[data];
data={{0,1}, {0.5, 2}, {1,0}, {1.5, -1}, {2, 1}}

N[DividedDifferece[data]];
TableForm[d]

Clear[P,x];
P[x_] = N[NewtonForm[data],10]

Clear[dots];
dots = ListPlot[data,PlotStyle->{PointSize[0.02]},
        DisplayFunction->Identity];
g = Plot[{P[x]}, {x, -0.2, 2.2},
DisplayFunction->Identity];
Show[g,dots, PlotLabel->"data and the Newton polynomial.",
AxesLabel->{"x","y"},DisplayFunction->$DisplayFunction];

Clear[f,x];
f[x_] = Cos[x]
xy = Table[{i,N[f[i],8]},{i,0,3}]

Clear[P3,x];
P3[x_] = N[NewtonForm[xy],10]

TableForm[d]

dots = ListPlot[xy,PlotStyle->{PointSize[0.02]},
        DisplayFunction->Identity];
pointset = Map[Point,xy];
graphset = Prepend[pointset,PointSize[0.02]];
dots = Graphics[graphset];
Print["The function is f[x] = ",f[x]];
g3 = Plot[ f[x], {x,-1, 4}, 
     DisplayFunction->Identity];
g4 = Plot[{P3[x]},{x,-1, 4},
PlotStyle -> Thickness[0.01],
DisplayFunction->Identity];
Show[g3,g4,dots,PlotRange->{{-1, 4},{-1.1,1.1}},
PlotLabel->"f(x) and the Newton polynomial.",
AxesLabel->{"x","y"},DisplayFunction->$DisplayFunction];