A continuous function on $[0,1]$ orthogonal to each monomial of the form $x^{n^2}$ Let us consider the continuous functions over $[0,1]$ fulfilling
$$ \int_{0}^{1} f(x) x^n\,dx = 0 $$
for $n=0$ and for every $n\in E\subseteq\mathbb{N}^+$. The Müntz–Szász theorem gives that
$$ \sum_{n\in E}\frac{1}{n} = +\infty \Longleftrightarrow f(x)\equiv 0 $$
so there is a non-zero continuous function $f(x)$ such that
$$ \int_{0}^{1} f(x) x^{n^2}\,dx = 0 \tag{1}$$
holds for every $n\in\mathbb{N}$.

Question: can we construct a nice, explicit function $f\neq 0$ fulfilling $(1)$ for every $n\in\mathbb{N}$?

We may consider functions of the form
$$ f(x) = \sum_{n\geq 0} c_n P_n(2x-1) $$
with $P_n(2x-1)$ being the $n$-th shifted Legendre polynomial. The orthogonality to $1$ and $x$ translates into $c_0=c_1=0$, the orthogonality to $x^4$ translates into $\frac{2}{35}c_2+\frac{1}{70}c_3+\frac{1}{630}c_4=0$, the orthogonality to $x^9$ translates into $\frac{3}{55}c_2+\frac{21}{715}c_3+ \frac{9}{715}c_4 +\frac{3}{715}c_5+\frac{3}{2860}c_6+\frac{9}{48620}c_7+\frac{1}{48620}c_8+\frac{1}{923780}c_9=0$ and so on. The minimal (with respect to the $\ell^2$ norm) solution of this infinite system with $c_2=1$ (or $c_4=1$) should give a sequence $\{c_n\}_{n\geq 0}$ ensuring the continuity of $f(x)$, but this is non-trivial and I would appreciate a more explicit construction / example of such $f$.
Addendum: another possible construction is to apply the Gram-Schmidt process to $\{1,x,x^4,x^9,\ldots\}$ in order to get a sequence of polynomials $\{p_n(x)\}_{n\geq 0}$ such that

*

*$p_n(x)=\sum_{k=0}^{n} c_k x^{k^2}$

*$n\neq m \Longrightarrow \langle p_n(x),p_m(x)\rangle = 0$

*$\max_{x\in [0,1]} |p_n(x)| = 1$ or $\langle p_n(x),p_n(x)\rangle = 1$
then take $f(x)$ as the pointwise limit of a convergent subsequence of $\{p_n(x)\}_{n\geq 0}$. Still, not really explicit.
A more promising approach is to consider some lacunary Fourier series, like
$$ g(\theta) = \sum_{n\geq 1}\frac{\cos(n\theta)}{n^2} - \sum_{n\geq 1}\frac{\cos(n^2\theta)}{n^4}, $$
which clearly fulfills $\int_{-\pi}^{\pi}g(\theta)\cos(n^2\theta)\,d\theta = 0$, then turn such $g(\theta)$ into an $f(x)$ fulfilling $(1)$ via some slick substitution.
Yet another way is to consider the inverse Laplace transform of
$$ \frac{1}{s}\prod_{k=0}^{n}\frac{k^2+1-s}{k^2+1+s} $$
evaluated at $-\log x$. This gives a polynomial, bounded between $-1$ and $1$, which is orthogonal to $1,x,x^4,\ldots,x^{n^2}$. Is this sequence of polynomials (or a subsequence of this sequence) convergent to a continuous function? I do not know. If so,
$$ f(x)=\mathcal{L}^{-1}\left(\frac{\sin(\pi\sqrt{s-1})}{\sqrt{s-1}}\cdot\frac{\sqrt{s+1}}{\sin(\pi\sqrt{s+1})}\cdot \frac{1}{s}\right)(-\log x)$$
is an explicit solution. Here it is a plot of the first polynomials produced by the last approach:

 A: The M-S theorem does not say what you say it says. (In a comment I pointed out that there is an important quantifier missing from your version - the current problem is maybe more serious...)
If $\sum_{n\in S}1/n<\infty$ the theorem says that $V=span\{x^n:n\in S\}$ is not dense. That implies that there is a non-zero real measure "orthogonal" to everything in $V$, but not a continuous function.
A: HINT:
We can write the perpendicular vector drawn from $v$ to the subspace spanned by the linearly independent vectors $v_1$, $\ldots$, $v_n$ as
$$\frac{D}{ G(v_1, \ldots, v_n)}$$
where $G(v_i)$ is the Gram determinant of the system $v_i$, while
$$D = \left | \begin{matrix} v&v_1& \ldots & v_n \\ \langle v, v_1\rangle  & \langle v_1, v_1 \rangle & \ldots & \langle v_n, v_1 \rangle\\ \ldots & & \ldots \\ \langle v, v_n\rangle & \langle v_1, v_n \rangle & \ldots & \langle v_n, v_n \rangle \end{matrix} \right|$$
Indeed, it is easy to see that $\langle D, v_i\rangle = 0$ for all $i$, and the coefficient of $v$ in $D$ equals $G(v_1, \ldots, v_n)$
Note also that
$$d(v, (v_1, \ldots, v_n))^2 = \frac{G(v, v_1, \ldots, v_n)}{G(v_1, \ldots, v_n)}$$
Now consider in $L^{2}[0,1]$ the vectors $v_i$ being the function $x\mapsto x^{\alpha_i}$, while $v$ is the function $x\mapsto x^{\alpha}$. We get to calculate some minors that turn out to be Cauchy determinants. Moreover, if $\sum\frac{1}{\alpha_i} < \infty$, and $\alpha \not \in \{\alpha_i\}_i$, then
$$d^2(x^{\alpha}, (x^{\alpha_i})_i) >0$$
being expressed as an (infinite) convergent product. We also get the perpendicular from $x^{\alpha}$ to $(x^{\alpha_i})_i$ as a limit of a  sequence convergent in $L^2[0,1]$.
$\bf{Added:}$ It's worth writing down the formulas for the projection of a vector $v$ onto the span of linearly independent $v_1$, $\ldots$, $v_n$
$$\pi_{\langle v_1, \ldots, v_n\rangle}(v) = \sum_{i=1}^n  \frac{ G( (v_1, \ldots, \overset{i}{v},
\ldots, v_n), (v_1, \ldots, v_n)) }{G((v_1, \ldots, v_n), (v_1, \ldots, v_n))}\cdot v_i $$
where by $G(u, v)$ we denot the (determinant of ) a matrix with entries $\langle u_k, v_l\rangle $.
Now, the Gram determinant $G((x^{\alpha_i}), (x^{\beta_j}))$ is the Cauchy determinant
$$C((\alpha_i), (\beta_j))\colon = \left|\frac{1}{\alpha_i + \beta_j+1}\right|$$ which equals
$$V(\alpha_i ) \cdot V(\beta_i) \cdot \prod_{i,j} \frac{1}{\alpha_i + \beta_j + 1}$$
where $V$ are the Vandermonde determinants.  We get the projection of the function $x^{\alpha}$ onto the span of $x^{\alpha_i}$
$$\sum_{i} \prod_{j, j\ne i} \frac{\alpha_j - \alpha}{\alpha_j - \alpha_i} \cdot \prod_j\frac{\alpha_j + \alpha_i + 1}{\alpha_j + \alpha+ 1} \cdot x^{\alpha_i}$$
Let's assume that $\alpha_j \ne 0$ for all $j$. We get the expression
$$\sum_i \frac{\prod_{j, j\ne i} (1- \frac{\alpha}{\alpha_j})}{\prod_{j, j\ne i} (1- \frac{\alpha_i}{\alpha_j})}\cdot \frac{\prod_{j} (1+ \frac{\alpha_i+1}{\alpha_j})}{\prod_{j} (1+ \frac{\alpha+1}{\alpha_j})}\cdot x^{\alpha_i}$$
Forecasting the case $\alpha_n = n^2$, let's consider the function
$$\phi(t) = \prod_j (1- \frac{t}{\alpha_j})$$
We have the projector formula
$$\sum_i \frac{\phi(\alpha)}{(1-\frac{\alpha}{\alpha_i}) \cdot (- \alpha_i \cdot \phi'(\alpha_i))} \cdot \frac{\phi(-(\alpha_i + 1))}{\phi(-(\alpha+1))} \cdot x^{\alpha_i}$$
Consider now the case $\alpha_k = k^2$. Then we get
$$\phi(t) = \prod_{k\ge 1} (1- \frac{t}{k^2}) = \frac{\sin \pi \sqrt{t}}{\pi \sqrt{t}}$$
while
$$\phi(-t) = \frac{\sinh \pi \sqrt{t}}{\pi \sqrt{t}}$$
One checks that at every $\alpha_k = k^2$ we have $ -k^2 \phi'(k^2) = \frac{(-1)^{k+1}}{2}$.
Now one puts all of it together. We get an analytic function.  More details later.
$\bf{Added:}$ The projector formula can also be written as
$$\sum_i \frac{\phi(\alpha)/(\alpha - \alpha_i)}{\phi'(\alpha_i)} \cdot \frac{\phi( -( \alpha_i + 1))}{\phi(-(\alpha+1))} \cdot x^{\alpha_i}$$
