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).
I managed to figure this out a little late.
In the given functional equation, replace $x$ with $-x$.
The equation is now written as:
$f(2x^2-1)=-2xf(-x)$
But it is given that $f(2x^2-1)=2xf(x)$
Therefore, $-2xf(-x)=2xf(x)$
$=>$ $-f(-x)=f(x)$
A rather simple functional equation which can be easily differentiated 4 times.
The answer is $0$.
Best Answer
I think the OP's question may really have been how to incorporate all the $f^{(n)}\text{'s}$ into a single formula, since $f, f',$ and $f''$ don't seem to follow the same pattern as the later ones.
Usually the answer would simply be written in cases, like
$$f^{(n)}(x)=\begin{cases}x^2 \log x,&\text{if }n=0, \\2x \log x+x, &\text{if }n=1, \\2\log x +3, &\text{if }n=2, \\2(-1)^{n+1}(n-3)!x^{2-n},&\text{if }n\ge 3, \end{cases}$$
without trying to come up with some artificial formula that covers $n=0, 1,$ and $2$ also.
If you really want to try to merge these, the best way would be to write
$$f^{(n)}(x)=\begin{cases} \frac{x^{2-n}}{(2-n)!}(2 \log x - 2 H_{2-n}+3),&\text{if }n\le 2, \\2(-1)^{n+1}(n-3)!x^{2-n},&\text{if }n\ge 3, \end{cases}$$
where $\,H_n=\sum_{k=1}^n \frac1{k}\,$ is the $\,n^{\text{th}}$ harmonic number, because then the formula even works for negative integer values of $n,$ with $\,f^{(-n)}(x)\,$ being an $n^{\text{th}}$ antiderivative of $\,f(x).$ (This still has two cases though, and I don't see a natural way to merge them into one.)
But if this was a homework problem, as I assume it was, this complicated formula isn't what was intended. And if you really only need the derivatives, it's overkill!