I am stuck on how to start. My idea is the following but I am not sure how it will lead to anything:
We know that because $f$ is continuous function on real interval $[a,b]$, then for all $\epsilon$, there is polynomial $p$ such that for all $x\in[a,b]$, we have that $\left | f(x)-p(x) \right |< \epsilon $ by Weierstrass' Theorem.
Construct Taylor series for $f(x)$, call it $p(x)$.
So $f(x)-p(x) =0 $ so it will clearly be less than any choice of $\epsilon$.
But, how do the conditions on $f$, that $f(0)=f(1)=0$ come into play? Am I even on right track?
Here it is a trick: since $f(x)$ is continuous on $[0,1]$, we have $|f(x)|\leq M$ for any $x\in[0,1]$.
Let $g(x)=\frac{1}{\sqrt{2M}}\sqrt{f(x)+2M}$, which is also a continuous function. By the Weierstrass approximation theorem there is some polynomial $g_n(x)$ such that $\left|g(x)-g_n(x)\right|\leq \frac{1}{n}$ for any $x\in[0,1]$. $$ h_n(x) = g_n(x)+(1-g_n(0))(1-x)+(1-g_n(1))x $$ is a polynomial such that $h_n(0)=h_n(1)=1$ and $|g(x)-h_n(x)|\leq \frac{2}{n}$ for any $x\in[0,1]$. $$ f_n(x) = 2M h_n(x)^2 - 2M $$ is a polynomial such that $f_n(0)=f_n(1)=0$ and $|f(x)-f_n(x)|\leq\frac{5M}{n}$ for any $x\in[0,1]$.
If $f(x)$ is differentiable at $0$ and $1$ we may simply approximate $\frac{f(x)}{x(1-x)}$ through polynomials $p_n(x)$, then consider $f_n(x)=x(1-x)p_n(x)$. Of course this approach does not work if $f(x)=\sqrt{x(1-x)}$, for instance. And in general not every continuous function agrees with its Taylor series at some point.