Poincare lemma for non-smooth differentiable forms The Poincare lemma is almost always formulated for differential forms with smooth coefficients (or sometimes for currents that have distributional coefficients). I would like to have it for $C^k$-coefficients. For one dimension everything is fine:
$$0\to \mathbb R \to C^1(\mathbb R) \to C(\mathbb R)\to 0$$
(where the first map assigns the constant function and the second is the derivative) is exact because of the fundamental theorem of calculus. For $\mathbb R^2$ the analogous sequence would be
$$0\to \mathbb R \to C^2(\mathbb R^2) \to C^1(\mathbb R^2)^2\to C(\mathbb R^2)\to 0$$
where the first map again assigns the constant function, the second is the gradient and the third is $d: C^1(\mathbb R^2)^2\to C(\mathbb R^2)$, $(g_1,g_2)\mapsto \partial_2 g_1 - \partial_1 g_2$.
This sequence is exact at the first three spots (the usual proof of the poincare lemma) but I don't see whether it is exact at the last spot, that is, whether $d$ is onto.
(For the same complex with $C^\infty$ (or $\mathscr D'$) instead of $C^k$ surjectivity of $d$ is trivial: For $h\in C^\infty(\mathbb R^2)$ one can just integrate one variable: $g_1(x,y)=\int_0^yh(x,t)dt$ and $g_2=0$ yield $d(g_1,g_2)=h$. However, if $h$ is just continuous there is no reason why $g_1$ should be differentiable w.r.t. $x$.)

EDIT. The discussion with Igor (and his updated answer) lead to the simple observation that it would be sufficient to show the surjectivity of $\Delta: C^2(\mathbb R^2) \to C(\mathbb R^2)$. I have therefore added the tags pde and fa. This last question looks so natural that the answer should be known (and I believe it to be negative).
 A: Roughly, any regularity/support class of forms that can be smoothed in a translation invariant way, without leaving the same regularity/support class, still satisfies the Poincaré lemma. An explicit example, which uses a homotopy operator similar but not identical to the one in your question, can be found in de Rham's original book (Ch.3, §15), where he discusses smoothing of currents (distributional forms) to $C^\infty$ forms.
Smoothing can be achieved by convolution with a $C^\infty$ compactly supported bump function. This smoothing map can then be show to be homotopic to the inclusion $C^\infty \subset C^k \subset \mathscr{D}'$. Thus, the de Rham cohomology in $C^k$ is isomorphic to the one in $\mathscr{D}'$, which in turn is isomorphic to the one in $C^\infty$. So the Poincaré lemma holds for all of them.

Update: This is just to write in a bit more detail what I mentioned in the comments about using the inverse Laplacian to construct a contracting homotopy. I'm essentially appealing to a variant of Hodge theory. Let $\Delta = \partial^i\partial_i$ denote the standard Laplacian on $\mathbb{R}^n$. I will denote the de Rham differential by $d$ and the co-differential by $\delta$, with explicit formulas $d[h]_{i_1\cdots i_k} = k\partial_{[i_1} h_{i_2\cdots i_k]}$ and $\delta[h]_{i_1\cdots i_k} = \partial^i h_{i i_1\cdots i_k}$. They satisfy the well-known identity $\Delta = \delta d + d\delta$, which means that $\delta$ can be thought of as a homotopy operator for the de Rham complex, inducing $\Delta$ as a null-homotopy morphism of the de Rham complex into itself.
Now, let $\Delta^{-1}$ be a Green function for $\Delta$. In other words, it is a translation invariant distribution $\Delta^{-1}(x-y)$ on $\mathbb{R}^n \times \mathbb{R}^n$ that is smooth away from the diagonal and satisfies $\Delta_x \Delta^{-1}(x-y) = \delta^n(x-y)$ (Dirac delta). We can use the formula $\Delta^{-1}[h]_{i_1\cdots i_k}(x) = \int \Delta^{-1}(x-y) h_{i_1\cdots i_k}(x)$ to make it act on forms with compact support. By translation invariance, both $\Delta$ and $\Delta^{-1}$ commute with $d$ and $\delta$.
To extend the action of $\Delta^{-1}$ on forms of non-compact support, we need a smooth compactly supported bump function $\phi$ that is identically $1$ on a neighborhood $U\ni 0$ in $\mathbb{R}^n$. Later we'll also need a partition of unity $\sum_a \chi_a = 1$ on $\mathbb{R}^n$, requiring that the supports of the $\chi_a$ are sufficiently small for each to fit into $U$ (perhaps with a safety margin) after a translation.
Let1 $\Delta^{-1}_\phi[h](x) = \int \phi(x-y) \Delta^{-1}(x-y) h(y) dy$. This operator now has the following important properties: (a) it is still translation invariant, so $d\Delta^{-1}_\phi = \Delta^{-1}_\phi d$; (b) $\Delta^{-1}_\phi[h]$ is well defined even if $h$ does not have compact support, since the defining integral now has compact support for any fixed $x$ because of the factor of $\phi(x-y)$; (c) $\Delta^{-1}_\phi[\Delta[h]] = h$ for any $h$ with sufficiently small support, which applies in particular to any argument of the form $h=\chi_a g$.
Finally, let us define the homotopy operator2 $H[h] = \Delta^{-1}_\phi[\delta h]$. It is now straight forward to verify the following identity for any $h$ with $d h = 0$:
$$
  d H[h] = \Delta^{-1}_\phi[d\delta h + \delta (d h)]
    = \Delta^{-1}_\phi[\Delta[h]]
    = \sum_a \Delta^{-1}_\phi[\Delta[\chi_a h]]
    = \sum_a \chi_a h
    = h .
$$
I think this actually answers the question, provided that $\Delta^{-1}[C^k_0] \subseteq C^{k+2}$, because then $H[C^k] \subset C^{k+1}$ and the homotopy works on the complex in question. Unfortunately, I'm not sure whether that's true or not. I think it is true provided one uses some $C^{k,\alpha}$ Hölder class instead of $C^k$, though.

1 Key phrase: properly supported parametrix for the Laplacian. 
2 Key phrase: properly supported parametrix for the de Rham complex.
A: A proof of the Poincaré lemma with optimal regularity for (non-integer order) Hölder and (nonnegative-order $L^p$, $2\leq p<\infty$) Sobolev forms is provided by Theorem 8.3, pp. 148-149 of the book by G. Csató, B. Dacorogna and O. Kneuss, "The Pullback Equation for Differential Forms" (Birkhäuser, 2011). It states in the Hölder case ($k\geq 0$ an integer, $0<\alpha<1$) that the sequence $$ 0\rightarrow\mathscr{C}^{k+n,\alpha}(U,\wedge^0(\mathbb{R}^n))\stackrel{d}{\rightarrow}\mathscr{C}^{k+n-1,\alpha}(U,\wedge^1(\mathbb{R}^n))\stackrel{d}{\rightarrow}\cdots\stackrel{d}{\rightarrow}\mathscr{C}^{k+1,\alpha}(U,\wedge^{n-1}(\mathbb{R}^n))\stackrel{d}{\rightarrow}\mathscr{C}^{k,\alpha}(U,\wedge^n(\mathbb{R}^n))\rightarrow 0 $$ is exact for any $U\subset\mathbb{R}^n$ open, bounded and contractible with smooth boundary (say, an open ball). The argument is actually quite similar to the one outlined by Igor in his answer (i.e. it uses a sharp version of the Hodge decomposition theorem), with appropriate care for the regularity gain in the above function spaces. It seems that it is also possible to provide a proof along more classical lines - the book quotes the PhD thesis of the first author for that, but I haven't got access to the thesis to check it yet.
The sequence made by the first two $d$ arrows even happens to be exact for $\alpha=0,1$ with the convention that $\mathscr{C}^{k,0}(U,\wedge^p(\mathbb{R}^n))=\mathscr{C}^k_b(\overline{U},\wedge^p(\mathbb{R}^n))$, by invoking instead Corollary 8.6, pp. 149-150 of the same book. One can derive from the above the exactness of the sequence $$ 0\rightarrow\mathscr{C}^{k+n,\alpha}_{loc}(\mathbb{R}^n,\wedge^0(\mathbb{R}^n))\stackrel{d}{\rightarrow}\mathscr{C}^{k+n-1,\alpha}_{loc}(\mathbb{R}^n,\wedge^1(\mathbb{R}^n))\stackrel{d}{\rightarrow}\cdots\stackrel{d}{\rightarrow}\mathscr{C}^{k+1,\alpha}_{loc}(\mathbb{R}^n,\wedge^{n-1}(\mathbb{R}^n))\stackrel{d}{\rightarrow}\mathscr{C}^{k,\alpha}_{loc}(\mathbb{R}^n,\wedge^n(\mathbb{R}^n))\rightarrow 0 $$ if one uses an exhaustion of $\mathbb{R}^n$ by (say) closed balls centered at the origin, since Theorem 8.3 guarantees that for all closed $\eta\in\mathscr{C}^{k,\alpha}(U,\wedge^r(\mathbb{R}^n))$ there is $\omega\in\mathscr{C}^{k+1,\alpha}(U,\wedge^{r-1}(\mathbb{R}^n))$ and a constant $C_{k,\alpha,U}>0$ such that $$d\omega=\eta\quad\text{with}\quad\|\omega\|_{\mathscr{C}^{k+1,\alpha}(U,\wedge^{r-1}(\mathbb{R}^n))}\leq C_{k,\alpha,U}\|\eta\|_{\mathscr{C}^{k,\alpha}(U,\wedge^{r}(\mathbb{R}^n))}$$ for all $r=1,\ldots,n$, $U$ as above. One can even have $\eta,\omega$ as above with $\alpha=0,1$ if $r=1$, so exactness of the first two $d$ arrows in the second sentence also holds with $\alpha=0,1$, despite the negative result regarding existence of $\mathscr{C}^2$ solutions of $\Delta u=f$ with continuous $f$ in my comment. The book does not mention the possibility of getting exactness for the whole first sequence with $\alpha=0,1$, but an inspection of the proof of Corollary 8.6 shows that one only needs Theorem 8.3 in the $L^p$ Sobolev case with $p>n$ sufficiently large together with Morrey's inequality in order to upgrade weak derivatives to classical derivatives. Therefore, I don't really see an obstacle for having exactness of the whole sequences above with $\alpha=0$ - if I understand correctly the argument, one just needs an appropriate version of Morrey's inequality for vector-valued functions instead.
A: The answer to the question is in fact NO. It is quite surprising for me that it is relatively recent. Theorem 5 in the article Additional regularity for solutions of PDE
by David Preiss (J. reine angew. Math. 485 (1997), 197—207) states that there is a continuous function $\psi:\mathbb R^2 \to\mathbb R$ with compact support such
that the equation
div$\Psi(x) = \psi (x) $ for almost every $x\in\mathbb R^2$
possesses no locally Lipschitz solution $\Psi:\mathbb R^2\to\mathbb R^2$.
I thank Pedro Lauridsen Ribeiro and Igor Khavkine for their very informative comments. In particular, the book mentioned by Pedro eventually lead me to the article of Preiss.
(Unfortunately, I cannot split the bounty.)
