Is the image of a Borel subset of $[0,1]$ under a differentiable map still a Borel set?
Even for $C^\infty$ functions $f\colon[0,1]\to[0,1]$, it is possible that $f(E)$ is not Borel for some Borel sets $E$.
Starting with the simpler case where $f$ is only required to be continuous, one way to construct a counterexample is to start with a space filling curve $g\colon[0,1]\to[0,1]^2$ whose image is the entire unit square. Then write $f=\pi\circ g$ where $\pi\colon[0,1]^2\to[0,1]$ is the standard projection $\pi(x,y)=x$. It is then a standard fact that there exist Borel sets $F\subset[0,1]^2$ whose projection $\pi(F)$ is not Borel. In fact, the projections of Borel sets are precisely the analytic sets, and there do exist analytic sets which are not Borel. Then, taking $E=g^{-1}(F)$, this is Borel set whose image under $f$ is the non-Borel set $\pi(F)$.
Now, moving on to the differentiable - and, even, smooth - case, the first thing to note is that the image of a differentiable function $g\colon[0,1]\to[0,1]^2$ must have zero Lebesgue measure, so cannot fill the unit square. However, this is not necessary. All we need is that the image of $g$ contains $C^2$ for some closed uncountable set $C\subset[0,1]$ (e.g., a Cantor set). Such sets are always Borel isomorphic to the interval $[0,1]$ (as are all uncountable Polish spaces). So, we can procede as above and let $E=g^{-1}(F)$ where $F\subset C^2$ is Borel such that $\pi(F)$ is not Borel.
So, it remains to find a closed uncountable $C\subset[0,1]$ and smooth $g\colon[0,1]\to[0,1]^2$ containing $C^2$ in its image. I will choose $C$ to be a 'thin Cantor set'. By this, I mean that $C$ is constructed in the same way as the fat Cantor set, except we remove larger and larger proportions of the set at each stage of the construction instead of smaller and smaller proportions. Choose a sequence $\alpha_0,\alpha_1,\ldots\in(0,1/2)$ with $\alpha_n\to0$ as $n\to\infty$. Then define a sequence of closed sets as follows. Let $C_0=[0,1]$. For each $n\ge0$, once $C_n$ has been chosen (and is a finite union of disjoint closed intervals), we remove the centre $1-2\alpha_n$ proportion of each interval in $C_n$ to obtain $C_{n+1}$. Then, $C=\bigcap_nC_n$. Equivalently, $C$ is the set of points which can be written as $$ (1-\alpha_0)x_0+\alpha_0(1-\alpha_1)x_1+\alpha_0\alpha_1(1-\alpha_2)x_2+\alpha_0\alpha_1\alpha_2(1-\alpha_3)x_3+\cdots $$ for $x_0,x_1,\ldots\in\lbrace0,1\rbrace$. Note that if we had $\alpha_0=\alpha_1=\cdots=1/3$ then this is the usual Cantor middle thirds set.
Now, I'll describe the function $g$. This will be done in an inductive fashion using self-similarity. We can set $g=g_0$ where, for each $n\ge0$, $g_n\colon[0,1]\to[0,1]^2$ contains scaled and maybe reflected copies of another such map $g_{n+1}$ in each of the four corner sections $[0,\alpha_n]^2$, $[0,\alpha_n]\times[1-\alpha_n,1]$, $[1-\alpha_n,1]^2$, $[1-\alpha_n,1]\times[0,\alpha_n]$ of the unit square, and follows straight lines in between these regions. This is as in the image below.
Let $\theta\colon[0,1]\to[0,1]$ be a smooth function with $\theta(0)=0$, $\theta(1)=1$ and such that its derivatives all vanish at the endpoints. For example, $\theta(x)$ could be taken to be the integral of $\exp(-1/(x(1-x)))$ appropriately scaled. Then, to be precise, we define $g=g_0$ where $g_n=(g_{n,1},g_{n,2})\colon[0,1]\to[0,1]^2$ are defined inductively by $$ g_n(x)=\begin{cases} \alpha_n(g_{n+1,2}(7x),g_{n+1,1}(7x))&{\rm if\ } 0\le x\le 1/7,\\[.3em] (0,\alpha_n+(1-2\alpha_n)\theta(7x-1))&{\rm if\ } 1/7 < x < 2/7,\\[.3em] \alpha_n(g_{n+1,1}(7x-2),g_{n+1,2}(7x-2)-1)+(0,1)&{\rm if\ } 2/7\le x\le 3/7,\\[.3em] (\alpha_n+(1-2\alpha_n)\theta(7x-3),1-\alpha_n)&{\rm if\ } 3/7 < x < 4/7,\\[.3em] \alpha_n(g_{n+1,1}(7x-4)-1,g_{n+1,2}(7x-4)-1)+(1,1)&{\rm if\ } 4/7 \le x \le 5/7,\\[.3em] (1,1-\alpha_n-(1-2\alpha_n)\theta(7x-5))&{\rm if\ } 5/7 < x < 6/7,\\[.3em] \alpha_n(-g_{n+1,2}(7x-6),1-g_{n+1,1}(7x-6))+(1,0)&{\rm if\ } 6/7 \le x \le 1. \end{cases} $$ This does uniquely define $g$. Note that if $\tilde g_n$ was any other sequence of functions to the unit square satisfying the same recurrence equations, then we have $\lVert\tilde g_n - g_n\rVert=\alpha_n\lVert\tilde g_{n+1} - g_{n+1}\rVert$, where $\lVert\cdot\rVert$ is the supremum norm. Also, as they both map into the unit square, $\lVert \tilde g_n-g_n\rVert\le\sqrt{2}$. Putting these together, $$ \lVert\tilde g_n-g_n\rVert\le\alpha_n\alpha_{n+1}\cdots\alpha_{n+r-1}\sqrt{2} $$ for all $r\ge0$. Letting $r$ go to infinity shows that $\tilde g_n=g_n$.
We can also construct a sequence of smooth approximations to $g$. For each $N\ge0$ define a map $g_N^{(N)}(x)=(\theta(x),0)$. Then, inductively define $g_{N-1}^{(N)},g_{N-2}^{(N)},\ldots,g_0^{(N)}$ using the recurrence equation above. Note that, at each step, $g^{(N)}_n$ is a smooth function with $g^{(N)}_n(0)=(0,0)$ and $g^{(N)}_n(1)=(1,0)$, and whose derivatives vanish to all orders at the end points. Next, using $D^rg$ to represent the r'th order derivatives of a function $g$, we have $$ \lVert D^rg_n^{(N)}\rVert=7^r\max\left(\alpha_n\lVert D^rg_{n+1}^{(N)}\rVert,(1-2\alpha_n)\lVert D^r\theta\rVert\right) $$ for all $r\ge1$. From this, we can see that $\lVert D^rg_n^{(N)}\rVert$ is bounded by $(7^r\alpha_n)\cdots(7^r\alpha_{n+k-1})7^r\lVert D^r\theta\rVert$ for some $k\ge0$. As we chose $\alpha_n$ to be tending to zero, this is bounded by some constant $L_r$ independently of $n,k$. Next, for any natural numbers $n < M\le N$, the recurrence equation gives $$ \begin{align} \lVert D^rg_n^{(M)}-D^rg_n^{(N)}\rVert &= 7^r\alpha_n\lVert D^rg_{n+1}^{(M)}-D^rg_{n+1}^{(N)}\rVert\cr &=(7^r\alpha_n)(7^r\alpha_{n+1})\cdots(7^r\alpha_{M-1})\lVert D^rg_M^{(M)}-D^rg_M^{(N)}\rVert\cr &\le(7^r\alpha_n)(7^r\alpha_{n+1})\cdots(7^r\alpha_{M-1})2L_r \end{align} $$ Again since we chose $\alpha_n$ tending to zero, this shows that $\lVert D^rg_n^{(M)}-D^rg_n^{(N)}\rVert$ tends to zero as $M,N$ go to infinity. Therefore, the sequence of functions $g_n^{(N)}$ converges uniformly as $N\to\infty$, along with its derivatives to all orders. The limit, $g_n$, satisfies the defining recurrence equations above. So, $g$ is infinitely differentiable.
Note, finally, that by choosing $x=x_0/7+x_1/7^2+x_2/7^3+\cdots$ for appropriate $x_n\in\lbrace0,2,4,6\rbrace$ we can find $x$ so that $g(x)$ is equal to any desired element of $C^2$.