Many years back, I wrote this blog post describing Matlab/GNU Octave code to compute a Lagrange interpolation. However, a one liner is possible using Vandermonde matrices.

function coeff = lagrange(xpts, fpts)

  %% Lagrange interpolation fitting polynomial p(x) to function f(x)
  %%
  %% p(x) = a_n x^n + a_{n-1} x^{n-1} + a_{n-2} x^{n-2} + ... + a_1 x + a_0
  %%
  %% xpts -- column vector of sampling points
  %% fpts -- column vector with value of the function at the sampling points
  %% coeff -- column vector of polynomial coefficients in the form [a_n, a_{n-1}, a_{n-2}, ... , a_1, a_0]'

  coeff = vander(xpts) \ fpts;

This use of a Vandermonde matrix follows from the fact that the Lagrange interpolation problem is just polynomial fitting a set of data points. This means that the obtained polynomial coefficients are just a solution to a linear system of equations.

While this code is an elegant one-liner, note that solving the Vandermonde system is a problem of high condition number, and errors will quickly become very large for any non-trivial number of points.