On the reals $\mathbb{R}$, the maximum entropy distribution with a given mean and variance is the Gaussian distribution.
Let $\mu, \sigma > 0$. What is the maximum entropy distribution on the natural numbers $\mathbb{N} = \{0, 1 , 2, \ldots \}$ with mean $\mu$ and variance $\sigma^2$?
This follows (modulo any minor technical details I haven't checked) from the theory of exponential families. The main result there says that the distribution which maximizes entropy subject to constraints on moments lies in the exponential family with sufficient statistics corresponding to these moments.
More concretely, let $n$ be distributed according to the distribution $\pi$ on $\{0,1,\ldots\}$. You constrain $\pi$ to satisfy $\mathbb{E}_\pi n = \mu$ and $\mathbb{E}_\pi n^2 = \sigma^2 + \mu^2$ (using the standard relation between the variance and uncentered second moment). Then $\pi$ must be of the form \[ \pi(n) = \frac{1}{Z}\exp\left(\alpha n + \beta n^2\right)\text{ for all }n, \] where $\alpha,\beta\in\mathbb{R}$ and $Z = \sum_n \exp\left(\alpha n + \beta n^2\right)$ is the constant that ensures this distribution normalizes. (Actually we can say $\beta\leq 0$ since otherwise we'd have $Z=\infty$ and the distribution wouldn't be normalizable. Similarly $\alpha<0$ if $\beta = 0$.)
So the problem reduces to one of finding $\alpha$ and $\beta$ given $\mu$ and $\sigma^2$, i.e. we've gone from needing to find infinitely many values $\pi(n)$ to only two. In this case I have a feeling that there is no closed form for $Z$ and so finding an explicit expression for $\alpha$ and $\beta$ in terms of $\mu$ and $\sigma^2$ is unlikely. However, there are optimization techniques for "moment matching" which will let you approximate these numerically.
In the simpler case when only the mean is constrained, things work out more nicely. If you go through the same sort of procedure but without the $\beta n^2$ term, you'll get a $Z$ which you can sum explicitly and which is finite when $\alpha < 0$. You can then relate $\mu$ and $\alpha$. The result is a geometric distribution with parameter $\frac{1}{\mu}$.
A straightforward calculation with Lagrange Multipliers gives $$p_n=e^{A+Bn+Cn^2}$$ where the three unknowns $A,B$ and $C$ are chosen to satisfy the three conditions $\sum p_n=1$, $\sum np_n=\mu$, and $\sum n^2p_n-\mu^2=\sigma^2$ .
