Consider the expression $$ f(x,n)=\sum_{k=0}^n\frac{e^{-k^2x}-e^{-(k+1)^2x}}{2k+1} $$ What is $\lim_{n\to\infty}f(x,n)$? If not a closed formula, is it possible to find an upper bound?
Note: This seems like some type of telescoping series, weighted by terms of the form $(2k+1)^{-1}$. The main reason I ask is that, if plotting $f$ for increasing values of $n$, we get what seems like an approximation to a limiting function of $x$, as seen here
Any ideas?

$$f(x)=\sum_{n=0}^\infty\frac{e^{-n^2x}-e^{-(n+1)^2x}}{2n+1}=\sum_{n\in\mathbb{Z}}\frac{e^{-n^2 x}}{2n+1}=\sum_{n\in\mathbb{Z}}\frac{e^{-n^2 x}}{1-4n^2}$$ is clearly related to the family of theta functions; we may express it in terms of $$\vartheta(x)=\sum_{n\in\mathbb{Z}}e^{-\pi(nx)^2},$$ chosen for its well-known property $\vartheta(1/x)=x\vartheta(x)$.
The last sum for $f(x)$ above (holds for $x\geqslant 0$ and) gives $$f(x)+4f'(x)=\sum_{n\in\mathbb{Z}}e^{-n^2 x}=\vartheta(\textstyle\sqrt{x/\pi})\\\implies f(x)=e^{-x/4}\int_0^{x/4}e^y\,\vartheta(2\textstyle\sqrt{y/\pi})\,dy.$$
In particular, for small $x$ we have $f(x)=\sqrt\pi D_+(\sqrt{x}/2)+O(xe^{-\pi^2/x})$, where $$D_+(z)=e^{-z^2}\int_0^z e^{t^2}\,dt=\frac12\sum_{n=0}^\infty\frac{(-1)^n n!}{(2n+1)!}(2z)^{2n+1}.$$ A less accurate estimate is then $f(x)=\sqrt{\pi x}/2+O(x^{3/2})$.