A closed subspace of $C([0,1])$ with all functions of bounded variation has finite dimension

Using somewhat advanced results from the theory of Banach spaces (in particular, Dvoretzky's theorem), one can make the idea suggested by @Gio67 for constructing functions $h \in E$ with total variation much bigger than their $L^\infty$-norm into a rigorous proof.

Let us assume towards a contradiction that there exists an infinite-dimensional, closed subspace $E \subset C([0,1])$ satisfying $E \subset BV$. As noted by the OP and also shown in the answer by @Gio67, there then exists a constant $M > 0$ satisfying $\mathrm{Var}(f) \leq M \cdot \| f \|_\infty$ for all $f \in E$.

In the following, consider $E$ as a Banach space equipped with the norm $\| \cdot \|_\infty$. Then, Dvoretzky's theorem (in the version given in this paper) shows for any given $N \in \mathbb{N}$ that there exists a subspace $V_N \subset E$ of dimension $\dim V_N = 2 N$ and such that $d(V_N, \ell_2^{2N}) < 2$, meaning that there exists an isomorphism $T : V_{N} \to \ell_2^{2N}$ satisfying $\| T \| \cdot \| T^{-1} \| < 2$. Here, $\ell_2^{2N}$ is $\mathbb{R}^{2N}$, equipped with the standard Euclidean norm.

Fix a basis $(f_1,\dots,f_{2N})$ for $V_N$. Then, Theorem 2 in this paper shows that there exist $x_1,\dots,x_{2N} \in [0,1]$ (so-called Fekete points) and functions $g_1,\dots,g_{2 N} : [0,1] \to \mathbb{R}$ with $\| g_i \|_{\infty} \leq 2$ and such that $$ f = \sum_{\ell=1}^{2 N} f(x_\ell) \, g_\ell \qquad \forall \, f \in V_{N}. \tag{$\ast$} $$

It is not too difficult to see that $$ g_1,\dots,g_{2N} \in V_N \quad \text{and} \quad g_i (x_k) = \delta_{i,k}. \tag{$\lozenge$} $$ Indeed, first note that $(\ast)$ implies $V_N \subset \mathrm{span} \{ g_1,\dots,g_{2N} \}$. Since $\dim V_N = 2 N$, this easily implies $V_N = \mathrm{span} \{ g_1,\dots,g_{2N} \}$ and in particular $g_1,\dots,g_{2N} \in V_N$. Hence, the linear maps $S : V_N \to \mathbb{R}^{2N}, f \mapsto (f(x_i))_{1 \leq i \leq 2N}$ and $R : \mathbb{R}^{2 N} \to V_N, c \mapsto \sum_{\ell=1}^{2N} c_\ell \, g_\ell$ are well-defined. By $(\ast)$, we see $R S = \mathrm{id}_{V_N}$. Therefore, $R$ is surjective and $S$ is injective. Since all dimensions agree, this means that $R,S$ are both bijective with $S = R^{-1}$ and hence $S R = \mathrm{id}_{\mathbb{R}^{2N}}$. This easily implies $g_i (x_k) = \delta_{i,k}$.

Equation $(\lozenge)$ implies in particular that $x_i \neq x_j$ for $i \neq j$. Hence, by relabeling, we can assume that $0 \leq x_1 < x_2 < \dots < x_{2N} \leq 1$. Now, let $(\epsilon_i)_{i=1,\dots,N}$ be independent Rademacher random variables (that is, $\mathbb{P}(\epsilon_i = 1) = \mathbb{P}(\epsilon_i = -1) = \frac{1}{2}$).

The idea is to use these Rademacher random variables and the fact that $V_N$ is "almost a Hilbert space" to introduce a suitable form of "almost orthogonality", in order to produce a function with $L^\infty$ norm of order $O(\sqrt{N})$ but with total variation of order $O(N)$.

Explicitly, note that if $h_1,\dots,h_N$ are elements of a Hilbert space $H$, then $$ \mathbb{E} \Big\| \sum_{i=1}^N \epsilon_i \, h_i \Big\|_H^2 = \sum_{i,j = 1}^N \Big( \langle h_i,h_j \rangle_H \cdot \mathbb{E} [\epsilon_i \epsilon_j] \Big) = \sum_{i=1}^N \| h_i \|_H^2 , $$ since $\mathbb{E} [\epsilon_i \epsilon_j] = \mathbb{E} [\epsilon_i] \cdot \mathbb{E} [\epsilon_j] = 0$ whenever $i \neq j$. Therefore, using the operator $T : V_N \to \ell_2^{2N}$ that we got from Dvoretzky's theorem, we see \begin{align*} \mathbb{E} \Big\| \sum_{i=1}^N \epsilon_i \, g_{2 i} \Big\|_\infty^2 & = \mathbb{E} \Big\| T^{-1} \sum_{i=1}^N \epsilon_i \, T g_{2 i} \Big\|_\infty^2 \\ & \leq \| T^{-1} \|^2 \cdot \mathbb{E} \Big\| \sum_{i=1}^N \epsilon_i \, T g_{2 i} \Big\|_{\ell_2^{2N}}^2 \\ & = \| T^{-1} \|^2 \cdot \sum_{i=1}^N \| T g_{2 i} \|_{\ell_2^{2N}}^2 \\ & \leq \| T^{-1} \|^2 \cdot \| T \|^2 \cdot \sum_{i=1}^N \| g_{2 i} \|_{\infty}^2 \leq 16 \cdot N , \end{align*} since $\| T^{-1} \| \cdot \| T \| < 2$ and $\| g_i \|_\infty \leq 2$. In particular, there exists at least one choice of $\epsilon_1,\dots,\epsilon_N \in \{ \pm 1 \}$ such that $$ h := \sum_{i=1}^N \epsilon_i \, g_{2 i} \in V_N $$ satisfies $$ \| h \|_\infty \leq \sqrt{16 \cdot N} = 4 \sqrt{N}, \quad \text{and hence} \quad \mathrm{Var}(h) \leq M \cdot \| h \|_\infty \leq 4 M \sqrt{N}. $$

To complete the proof, note that $(\lozenge)$ and the definition of $h$ imply that $h(x_\ell) = 0$ for $\ell$ odd and that $h(x_{2 i}) = \epsilon_i$ for $1 \leq i \leq N$. Therefore, $$ \mathrm{Var}(h) \geq \sum_{\ell=1}^{2N - 1} \bigl|h(x_{\ell+1}) - h(x_\ell)\bigr| \geq \sum_{i = 1}^{N} \bigl|h(x_{2 i}) - h(x_{2 i - 1})\bigr| = \sum_{i=1}^{N} |\epsilon_i| = N. $$ Hence, we see that $N \leq \mathrm{Var}(h) \leq 4 M \sqrt{N}$, where $M$ is fixed and $N \in \mathbb{N}$ can be chosen arbitrarily large. This is clearly impossible and thus provides the desired contradiction.

This is not a complete answer, just an idea for possible proof.

The space $BV([0,1])$ of functions of bounded variation is a Banach space with the norm $$ \Vert f\Vert_{BV}=|f(0)|+\operatorname*{Var}f. $$ Consider the space $Z=\{f\in C([0,1]):\,\operatorname*{Var}f<\infty\}$ endowed with the norm $$ \Vert f\Vert_{Z}=\Vert f\Vert_{\infty}+\operatorname*{Var}f. $$ It is a Banach space, since if we take a Cauchy sequence in $Z$, it is a Cauchy sequence in $BV([0,1])$ and in $C([0,1])$ and so it converges to a function in $Z$.

Now consider the linear function \begin{align*} T & :E\rightarrow Z\\ f & \mapsto f \end{align*} If $f_{n}\rightarrow f$ in $E$ and $T(f_{n})\rightarrow g$ in $Z$, then $f=g$. Hence, $T(f)=g$. It follows by the closed graph theorem that $T$ is continuous, that is, there exists $M>0$ such that $$ \Vert T(f)\Vert_{Z}\leq M\Vert f\Vert_{\infty}% $$ for all $f\in E$. In particular, $$ \operatorname*{Var}f\leq M\Vert f\Vert_{\infty}% $$ for all $f\in E$.

Now the idea would be to use the fact that $E$ has infinite dimension to find a function $f\in E$ such that $\Vert f\Vert_{\infty}=1$ and $\operatorname*{Var}% f>M$. If we consider $n$ points $x_{1},\ldots, x_{n}\in\lbrack0,1]$ we can find $f\in E$ such that $f\neq0$ and $f(x_{1})=\ldots=f(x_{n})=0$. Indeed, if not, we consider the function $L:E\rightarrow\mathbb{R}^{n}$ given by $L(f)=(f(x_{1}),\ldots,f(x_{n}))$. Then $L$ is linear and $\ker L=\{0\}$, which implies that $E$ has finite dimension.

Using this fact, we can construct by induction a sequence of functions $f_{n}$ and of points $x_{n}$, such that $\max|f_{n}|=|f_{n}(x_{n})|>0$ and $f_{n}(x_{k})=0$ for all $k=1,\dots,n-1$. By rescaling, we can assume that $\max|f_{n}|=1$. Now one should try to use this fact to construct a function in $E$ with $\Vert f\Vert_{\infty}=1$ and $n$ ordered points such that $f(x_{2k})=0$ and $f(x_{2k+1})=1$. I am not quite sure how.