f[x_]:=1/(1+x^2);
pts={-2,-1,0,1,2};

v=Transpose[Join[{ConstantArray[1,5]},Table[pts^i,{i,1,4}]]]
MatrixForm[v]

coeff=LinearSolve[v,f[pts]]

poly=Sum[coeff[[i]]*x^(i-1),{i,1,5}]

Plot[{f[x],poly},{x,-2,2}]

f[x_]:=1/(1+x^2);
a=-1;
b=2;
n=10;

pts=Table[a+i*(b-a)/n,{i,0,n}]

v=Transpose[Join[{ConstantArray[1,n+1]},Table[pts^i,{i,1,n}]]]
MatrixForm[v]

coeff=LinearSolve[v,f[pts]]

poly=Sum[coeff[[i]]*x^(i-1),{i,1,n+1}]

Plot[{f[x],poly},{x,a-1,b+1}]