Prove $|f(x)-p_{3}(x)| \leq \frac{1}{384} \max\limits_{0 \leq x \leq 1}|f^{(4)}(x)|$

First define $$ g(x) = f(x) - p_3(x) \\ q(x) = x^2(1-x)^2 $$ and note that $g(0) = g'(0) = g(1) = g'(1) = 0$ and $g^{(4)} = f^{(4)}$.

$q \not \equiv 0$ is chosen as a polynomial satisfying $q(0) = q'(0) = q(1) = q'(1) = 0$ of smallest possible degree.

Now let $0 < c < 1$ and define $$ h(x) = q(x) g(c) - q(c) g(x) \, . $$ Then $$ h(0) = h(c) = h(1) = 0 \\ h'(0) = h'(1) = 0 $$ and repeated application of Rolle's theorem shows that $h^{(4)}(\xi) = 0$ for some $\xi \in (0, 1)$. So $$ 0 = h^{(4)}(\xi) = q^{(4)}(\xi) g(c) - q(c) g^{(4)}(\xi) \\ = 24 (f(c)-p_3(c)) - c^2(1-c)^2 f^{(4)}(\xi) $$ and therefore $$ f(c)-p_3(c) = \frac{f^{(4)}(\xi)}{24} c^2(1-c)^2 \, . $$ Finally, using $c(1-c) \le 1/4$, we conclude that $$ |f(c)-p_3(c)| \le \frac{1}{24 \cdot 4 \cdot 4} \max_{0 \leq x \leq 1}|f^{(4)}(x)| \, . $$

Let $x_0 \ne 0,1$, and $$K = \frac{f(x_0) - p(x_0)}{x_0^2(x_0-1)^2}$$ Then the function $$r(x) = f(x) - p(x) - K x^2(x-1)^2$$ satisfies $$r(0) = r'(0) = 0\\ r(x_0) = 0\\ r(1) = r'(1) = 0$$ Hence $r$ has $5$ zeroes, counting multiplicities. Let's show that there exists $\xi \in (0,1)$ such that $r^{(4)}(\xi)=0$. This would be clear if the zeroes were distinct,but it also works here. Indeed, by Rolle, there exist $\alpha\in (0, x_0)$, $\beta\in (x_0, 1)$ such that $r'(\alpha) = r'(\beta) = 0$. Now we have $4$ distinct roots of $r'$ and we can finally get $\xi$ such that $r^{(4)}(\xi) =0$. But calculating we get $$r^4(x) = f^{(4)}(x) - 24 K$$ and so $K= \frac{f^4(\xi)}{24}$. We therefore get $f(x_0) - p(x_0) =\frac{f^{(4)}(\xi)}{24} x_0^2(x_0-1)^2$. Note that $\xi$ depends on $x_0$. Note that a similar argument would work for any $x_0$, not necessarily in $(0,1)$. Therefore, we have the error term formula for Hermite interpolation $$f(x) - p(x)= \frac{f^{(4)} (\xi_x)}{24}x^2 (x-1)^2$$ where $\xi_x$ lies between $x$ and $0$,$1$.

$\bf{Added:}$ The explicit form for the Lagrange interpolation polynomial and an alternate proof can be done using generalized Vandermonde determinant, for instance take a look at my answer here. You may obtain by passing to limit the remainder for Hermite interpolation. Or one can use confluent determinants, see for instance Turnbull- Theory of Equations. Notice the the quotient of determinants in the answer is the quotient of remainders for two functions. If the bottom function is $(x-x_1)\cdots(x-x_n)$ the estimate is for Lagrange interpolation. Now , the generalization is for quotients of remainders $$\frac{f(x) -p_f(x)}{g(x) - p_g(x)} = \frac{f^{(n)}(\xi_x)}{g^{(n)}(\xi_x)}$$ for interpolation at $n$ points. Note that if the points coincide, we have the Taylor polynomial approximation.