Is it possible to use polynomial interpolation to show that $\cos(x) = \sum_{n=0}^{\infty} \frac{(-1)^nx^{2n}}{(2n)!}$?

The least squares fit polynomial to a function $f(x)$ over an interval $[a,b]$ is the polynomial $p(x)=a_0 + a_1x + \dots + a_dx^d$ with coefficients $a_i$ chosen to minimize $$ \int_{a}^{b} (p(x) - f(x))^2 \, dx. $$

The least squares fit polynomial to $f(x)=\cos(x)$ over the interval $[-\pi,\pi]$ actually does approach the Taylor series polynomial as you increase the degree.

---

A faster and more numerically stable way to compute the least squares fit polynomial for $f(x)$ is to let $p(x)$ be a linear combination of orthogonal polynomials (such as the Legendre polynomials). In other words, let $p(x) = \sum_{i=0}^d a_i P_i(x),$ where $P_i(x)$ is the Legendre polynomial of degree $i.$ The solution for $a_i$ is then $$ a_i = \frac{2i+1}{2} \int_{-1}^1 P_i(x) f(x) \, dx. $$

You can do what you wanted to do here with Mathematica.

Unprotect[Power]; 0^0 = 1;
ClearAll[sys, doit, a];
sys[n_?OddQ] := Table[With[{x = Pi (m - (n - 1)/2)/(n - 1)}, 
   Sum[a[k] x^k/k!, {k, 0, n - 1}] == Cos[x]], {m, 0, n - 1}];
doit[n_?OddQ] := (soln = Solve[sys[n], Array[a, n, 0]][[1]];
   f[x_] := Sum[a[k] x^k/k!, {k, 0, n - 1}] /. soln);
doit[7]; Print[soln, N[Table[soln[[n, 2]], {n, 7}], 5]];
Plot[{Cos[x], f[x]}, {x, -Pi, Pi}]
(*
{a[0] -> 1, a[1] -> 0, a[2] -> (-517 + 270*Sqrt[3])/(5*Pi^2),
a[3] -> 0, a[4] -> (-216*(-68 + 39*Sqrt[3]))/Pi^4,
a[5] -> 0, a[6] -> (46656*(-26 + 15*Sqrt[3]))/Pi^6}
{1.0000, 0, -0.99996, 0, 0.99789, 0, -0.93361}
*)

The purpose of 0^0=1 is to avoid a glitch. The sys[] function returns a list of polynomial equations, The main doit[] function solves the polynomials equations to find the coefficients. The doit[7] sets the function f[] to be a 6th degree polynomial with the coefficients found by solving the system of equations. The found coefficients are then printed and the Plot[] does a comparison of the actual Cos[] function with the found f[] polynomial function.

The key idea here is to interpolate the $\,\cos\,$ function at a set of points uniformly distributed over the whole interval $\,[-\pi,\pi].$ The difference in the function and its polynomial interpolation is pretty good even for $\,n=7.$ In fact, the coefficients are getting very close to the actual coefficients of the Taylor series.