Derivatives of Lagrange polynomials It seems there is some relationship between Lagrange polynomial and Legendre polynomial. That is Lagrange polynomial can be expressed as a function of Legendre polynomial. If so, I could use this relationship to code up in Matlab. I have been searching the relation for a long time but I did not find it. 
My question is: Could anyone tell me the relationship between Lagrange and Legendre polynomial? 
I am trying to code up in Matlab to calculate any order derivative of Lagrange polynomial. I know how to calculate the first derivative like here . But this method need about N FOR loops to calculate Nth derivative.
Thank you all!

You can find coefficients of Lagrange interpolation polynomial or any of its derivatives relatively easy if you use a matrix form of Lagrange interpolation presented in "Beginner's guide to mapping simplexes affinely", section "Lagrange interpolation". General formula for the polynomial looks as follows
$$
f(x) = (-1)
\frac{
    \det
    \begin{pmatrix}
        0       & f_0       & f_1       & \cdots & f_n       \\
        x^n     & x_0^n     & x_1^n     & \cdots & x_n^n     \\
        x^{n-1} & x_0^{n-1} & x_1^{n-1} & \cdots & x_n^{n-1} \\
        \cdots  & \cdots    & \cdots    & \cdots & \cdots    \\
        x       & x_0       & x_1       & \cdots & x_n       \\
        1       & 1         & 1         & \cdots & 1         \\
    \end{pmatrix}
}{
    \det
    \begin{pmatrix}
        x_0^n     & x_1^n     & \cdots & x_n^n     \\
        x_0^{n-1} & x_1^{n-1} & \cdots & x_n^{n-1} \\
        \cdots    & \cdots    & \cdots & \cdots    \\
        x_0       & x_1       & \cdots & x_n       \\
        1         & 1         & \cdots & 1         \\
    \end{pmatrix}
}.
$$
Here $(x_0;f_0)$, $\dots$, $(x_n;f_n)$ are the points it passes through. Using Laplace expansion along the first column you can get expressions for coefficients at $x^i$.
DERIVATIVE
If we take derivative of the expression above, it will only act on the first column of the matrix in the numerator (the only one containing $x$'s). One can prove this by expanding determinant in the numerator along the first column, taking derivative and then collecting everything back. As a result, first derivative looks as follows
$$
f'(x) = (-1)
\frac{
    \det
    \begin{pmatrix}
        0             & f_0       & f_1       & \cdots & f_n       \\
        n x^{n-1}     & x_0^n     & x_1^n     & \cdots & x_n^n     \\
        (n-1) x^{n-2} & x_0^{n-1} & x_1^{n-1} & \cdots & x_n^{n-1} \\
        \cdots        & \cdots    & \cdots    & \cdots & \cdots    \\
        1             & x_0       & x_1       & \cdots & x_n       \\
        0             & 1         & 1         & \cdots & 1         \\
    \end{pmatrix}
}{
    \det
    \begin{pmatrix}
        x_0^n     & x_1^n     & \cdots & x_n^n     \\
        x_0^{n-1} & x_1^{n-1} & \cdots & x_n^{n-1} \\
        \cdots    & \cdots    & \cdots & \cdots    \\
        x_0       & x_1       & \cdots & x_n       \\
        1         & 1         & \cdots & 1         \\
    \end{pmatrix}
}.
$$
Higher-order derivatives can be calculated the same way --- you change the first column appropriately, all the rest remains.
CODE
I'm afraid, I don't know Matlab good enough to provide you with appropriate code, but I worked with Python a little bit, maybe the following can help (sorry for bad codestyle -- I'm mathematician, not programmer)
import numpy as np
# input
x = [0, 2, 4, 5]  # <- x's
y = [2, 5, 7, 7]  # <- y's
# calculating coefficients
M = [[_x**i*(-1)**(i*len(x)) for _x in x] for i in range(len(x))]
C = [np.linalg.det((M+[y]+M)[d:d+len(x)]) for d in range(len(x)+1)]
C = (C / C[0] * (-1)**(len(x)+1) )[1:]
# polynomial lambda-function
poly = lambda _x: sum([C[i] * _x**i for i in range(len(x))])
# output and tests
print("Coefficients:\n", C)
print("TESTING:")
for _x, _y in zip(x, y):
    result = "[OK]" if np.allclose(_y, poly(_x)) else "[ERROR]"
    print(_x, " mapped to: ", poly(_x), " ; expected: ", _y, result)

This code calculates coefficients of the Lagrange interpolation polynomial, prints them, and tests that given x's are mapped into expected y's. You can test this code with Google colab, so you don't have to install anything. Probably, you can translate it to Matlab and modify so that it computes derivatives (multiply every summand by appropriate constant and modify powers of $x$ accordingly).
