I really like the answer of Martin Argerami. The only point that worried me was
(this is not that easy to write or read, but it is very easy to check if you draw the squares)
To circumvent this problem, one can use the following argument that is based on Banach's fixed point theorem (note that the geometric idea is still the same as for the construction of the Hilbert curve and/or in the post of Martin Argerami). Also compare the following to the post of john mangual.
We define mappings
$$
\begin{eqnarray*}
A_{1}: & \left[0,1\right]^{2}\rightarrow\left[0,\frac{1}{2}\right]^{2}\subset\left[0,1\right]^{2}, & x\mapsto\frac{1}{2}\left[\left(\begin{array}{cc}
0 & 1\\
-1 & 0
\end{array}\right)\left(x-\left(\begin{matrix}1/2\\
1/2
\end{matrix}\right)\right)+\left(\begin{matrix}1/2\\
1/2
\end{matrix}\right)\right]\\
& & \phantom{x}=\frac{1}{2}\left(\begin{array}{cc}
0 & 1\\
-1 & 0
\end{array}\right)\,x+\left(\begin{matrix}0\\
1
\end{matrix}\right),\\
A_{2}: & \left[0,1\right]^{2}\rightarrow\left[0,\frac{1}{2}\right]\times\left[\frac{1}{2},1\right]\subset\left[0,1\right]^{2}, & x\mapsto\frac{1}{2}x+\left(\begin{matrix}0\\
1/2
\end{matrix}\right),\\
A_{3}: & \left[0,1\right]^{2}\rightarrow\left[\frac{1}{2},1\right]^{2}\subset\left[0,1\right]^{2}, & x\mapsto\frac{1}{2}x+\left(\begin{matrix}1/2\\
1/2
\end{matrix}\right),\\
A_{4}: & \left[0,1\right]^{2}\rightarrow\left[\frac{1}{2},1\right]\times\left[0,\frac{1}{2}\right]\subset\left[0,1\right]^{2}, & x\mapsto\frac{1}{2}\left[\left(\begin{array}{cc}
0 & -1\\
1 & 0
\end{array}\right)\left(x-\left(\begin{matrix}1/2\\
1/2
\end{matrix}\right)\right)+\left(\begin{matrix}1/2\\
1/2
\end{matrix}\right)\right]+\left(\begin{matrix}1/2\\
0
\end{matrix}\right)\\
& & \phantom{x}=\frac{1}{2}\left(\begin{array}{cc}
0 & -1\\
1 & 0
\end{array}\right)x+\left(\begin{matrix}1\\
0
\end{matrix}\right).
\end{eqnarray*}
$$
These are just the rotations (+scaling and shifting) that are used in the numbering as in the post of Martin Argerami.
It is easy to check that these are well-defined and surjective maps (even bijective) that satisfy
$$\tag {*}
\begin{equation}
\left\Vert A_{j}x-A_{j}y\right\Vert _{2}=\frac{1}{2}\cdot\left\Vert x-y\right\Vert _{2} \end{equation}
$$
for all $j=1,\dots,4$ and $x,y \in [0,1]^2$.
For $j=1,\dots,4$ we also set $I_{j}:=\left[\frac{j-1}{4},\frac{j}{4}\right]\subset\left[0,1\right]^{2}$ and
$$
\begin{eqnarray*}
\varphi_{1}: & I_{1}\rightarrow\left[0,1\right], & x\mapsto1-4x,\\
\varphi_{2}: & I_{2}\rightarrow\left[0,1\right], & x\mapsto4\cdot\left(x-\frac{1}{4}\right),\\
\varphi_{3}: & I_{3}\rightarrow\left[0,1\right], & x\mapsto4\cdot\left(x-\frac{1}{2}\right),\\
\varphi_{4}: & I_{4}\rightarrow\left[0,1\right], & x\mapsto1-4\cdot\left(x-\frac{3}{4}\right).
\end{eqnarray*}
$$
It is again an easy exercise to show that these are well-defined and surjective (even bijective).
Now, let
$$
M:=\left\{ f\in C\left(\left[0,1\right];\mathbb{R}^{2}\right) \,\middle|\,f\left(0\right)=\left(\begin{matrix}0\\
0
\end{matrix}\right),\; f\left(1\right)=\left(\begin{matrix}1\\
0
\end{matrix}\right)\;,f\left(\left[0,1\right]\right)\subset\left[0,1\right]^{2}\right\} .
$$
It is easy to see that this is a closed, nonempty subset of $C([0,1]; \Bbb{R}^2)$, when this space is equipped with the $\sup$-norm (let us use $|\cdot| = \Vert \cdot \Vert_2$ on $\Bbb{R}^2$).
We now define an operator $T : M \rightarrow M$ by
$$
\left(Tf\right)\left(x\right):=A_{j}\left(f\left(\varphi_{j}\left(x\right)\right)\right)\text{ for }x\in I_{j}\text{ and }j\in\{1,2,3,4\}
$$
One can check that $Tf$ is a well-defined, continuous function for $f \in M$ and even $Tf \in M$ for $f \in M$, so that $T$ is well-defined.
Using the estimate $(\ast)$ above, one easily gets
$$
\left\Vert Tf-Tg\right\Vert _{\sup}\leq\frac{1}{2}\cdot\left\Vert f-g\right\Vert _{\sup}
$$
for all $f,g \in M$. Thus, by Banach's fixed point theorem, there is a (unique) fixed point $f_0 \in M$ of $T$, which is thus a continuous map $f_0 : [0,1] \rightarrow [0,1]^2$. It remains to show that $f_0$ is surjective. As the image $f_0 ([0,1])$ is compact, it suffices to show that it is dense in $[0,1]^2$.
Using induction on $j \in \Bbb{N}_0$, we will show that for each $y \in [0,1]^2$, there is some $x \in [0,1]$ satisfying
$$\tag {**}
\left\Vert f_{0}\left(x\right)-y\right\Vert _{2}\leq\sqrt{2}\cdot2^{-j}.
$$
This is clear for $j=0$, as the diameter of $[0,1]^2$ is just $\sqrt{2}$.
For the induction step, note that there is some $j \in \{1, \dots, 4\}$ such that $y \in A_j ([0,1]^2)$, i.e. $y = A_j y'$ for some $y' \in [0,1]^2$. By induction, there is some $x' \in [0,1]$ such that $(\ast \ast)$ is satisfied for $x'$ and $y'$ instead of $x,y$.
Furthermore (because $\varphi_j$ is surjective), there is $x \in I_j$ with $x' = \varphi_j (x)$.
Now, the fixed point property $f_0 =Tf_0$ implies
$$
f_{0}\left(x\right)=\left(Tf_{0}\right)\left(x\right)=A_{j}\left(f_{0}\left(\varphi_{j}\left(x\right)\right)\right)=A_{j}\left(f_{0}\left(x'\right)\right),
$$
Using the estimate $(\ast)$, we arrive at
\begin{align}
\left\Vert f_{0}\left(x\right)-y\right\Vert _{2}&=\left\Vert A_{j}\left(f_{0}\left(x'\right)\right)-A_{j}y'\right\Vert _{2}=\frac{1}{2}\cdot\left\Vert f_{0}\left(x'\right)-y'\right\Vert _{2}\\
&\leq\frac{1}{2}\cdot\sqrt{2}\cdot2^{-j}=\sqrt{2}\cdot2^{-\left(j+1\right)},
\end{align}
which completes the proof.