(*Example III *) 

Clear[data];

Clear[fit];

F[x_]:=1/(1+25*x*x);(*Runge's function*)

(*range*)
a=-1.0;
b=1.0;

plotfunction=Plot[F[x],{x,a,b},PlotStyle->{Hue[1],Thickness[0.01]},PlotLabel->"Original function",DisplayFunction->Identity];

n=7;h=(b-a)/n;(*define sample points as values of F[x]*)

data=Table[N[{x,F[x]}],{x,a,b,h}]

(* data=Table[N[{0.5*((a+b)+(b-a)*Cos[(i/n)*Pi]),F[0.5*((a+b)+(b-a)*Cos[(i/n)*Pi])]}],{i,0,n}] *) 

plotdata=ListPlot[data,PlotLabel->"Data points",PlotStyle->{PointSize[Large],Black}];

Show[plotfunction,plotdata] 

fit=InterpolatingPolynomial[data,x]

plotfit=Plot[fit,{x,a,b},PlotStyle->{Hue[0.5],Thickness[0.005]},PlotLabel->"interpolation",DisplayFunction->Identity];

Show[plotfunction,plotdata,plotfit,PlotLabel->{"red - original f(x) ,
blue - interpolation, black - nodes (sample points) "},Prolog->AbsolutePointSize[5],DisplayFunction->$DisplayFunction, PlotRange -> All]


Print["------------------------------------------------------"]

Print[" Interpolating function, Pn = " Expand[Simplify[fit]]]

(*Possible fixes below*)
(*Just cut and paste into the appropriate place in the above example*)

(*-------------------------------------------------------*)

(*Fix 1. Increase n*)(*Generally BAD idea*)(*Fix 2. Extend range of data. Approximates well if deep inside max range*)(*Fix 3. Use non-uniform or "unequally spaced" data points ("nodes","sampling")*)(*e.g.n Chebyshev points xi=cos[Pi*(2i+1)/(2n+2)],i≤n*)(*data=Table[N[{0.5*((a+b)+(b-a)*Cos[(i/n)*Pi]),F[0.5*((a+b)+(b-a)*Cos[(i/n)*Pi])]}],{i,0,n}]*)
(*Fix 4. Use more appropriate MODEL (basis functions).Gaussians in OUR case*)(*fit=Fit[data,Table[Exp[-i*1.0*x^2],{i,1,n}],x]*)
(*Thanks to Niteesh Bharara CS3414 for fixing the multiple plot problem*)