(** DIFFERENZE DIVISE **)

Clear[f]
n=4;
x[0]=0;f[0]=-1.0;
x[1]=0.25;f[1]=0.5;
x[2]=0.5;f[2]=2.0;
x[3]=0.75;f[3]=3.0;
x[4]=1.0;f[4]=2.5;
(**x[0]=1.;f[0]=0.7651977;
x[1]=1.3;f[1]=0.6200860;
x[2]=1.6;f[2]=0.4554022;
x[3]=1.9;f[3]=0.2818186;
x[4]=2.2;f[4]=0.1103623;**)

ff=Table[{x[k],f[k]},{k,0,n}]

ticks1={{0,0.2,0.4,0.6,0.8,1},{-1,0,1.2,3}};
pl1=ListPlot[ff,PlotStyle->PointSize[0.02],Ticks->ticks1];

(* ALGORITMO 3.2 del Faires Burden *)
Do[d[k][0]=f[k],{k,0,n}]
Do[d[k][j]=(d[k][j-1]-d[k-1][j-1])/(x[k]-x[k-j]),{k,1,n},{j,1,k}]
pn[x_]=Sum[d[k][k]*Product[(x-x[j]),{j,0,k-1}],{k,0,n}];
pl2=Plot[pn[x],{x,0,1.2}];

Show[pl1,pl2,PlotRange->{{0,1.2},{-1,3.5}}]

Do[Print["k = ", k,"  x = ",x[k], "      f = ",f[k]],{k,0,n}]

Do[Print["k = ", k,"  j = ",j, "   d = ", d[k][j]],{k,0,n},{j,0,k}]