I don't have time to get all the details down, but maybe this answer can give some helpful intuition.
Basically, something that will work is taking the first $N$ binary binary of your number in $[0,1]$ and assigning them as the most significant digits for each of the $N$ dimensions (so, the first decimal digit is the most significant digit in the first dimension, the second decimal digit is the most significant digit in the second dimension, and so on). Then, you perform a change of coordinates in the $[0,1]^N$ interval you are mapping to involving translations/rotations/reflections that depends on which ones were the previously read $N$ digits, and repeat the procedure.
A procedure in the case of $N=2$ that looks to me like it is following this model is given here: http://bit-player.org/2013/mapping-the-hilbert-curve. The operators for the change of coordinates are what they call the $H_i$. The main computational difference with your situation is that instead of considering quadits (groups of two bits, or which is the same, digits from $0$ to $3$), you are considering groups of $N$ bits (or which is the same, digits from $0$ to $2^N - 1$. Your $H_i$ operators will operate on size $N$ vectors, not size $2$ ones. And there will be $2^N$ of the $H_i$ operators, so you will not be able to precompute all of them.
Hope this clarifies more than it confuses.
Also, how much of this is a problem depends on the function you are trying to plot, but an issue you are likely to face is that to get a single digit of precision in each of the N coordinates, you must have $N$ digits of precision in your real number between $0$ and $1$. And 10,000 digits of precision is quite a bit : )