Why two 'modes' for a Poisson Distribution with $\lambda$ an integer? Let $X\sim \operatorname{Poi}[\lambda]$ such that $\lambda\in\mathbb{N}$.
$$\mathbb{P}[X=\lambda-1]=\frac{e^{-\lambda}\lambda^{\lambda-1}}{(\lambda-1)!}=\frac{e^{-\lambda}\lambda^\lambda}{\lambda!}=\mathbb{P}[X=\lambda].$$
Can this be explained on an intuitive level?
 A: Although the result is easy to verify when one is given a formula for a Poisson distribution, that provides little for the intuition to work with.  Instead, what we need is to develop an understanding of this distribution in terms of some probabilistic mechanism that we can see, experiment with, and understand easily: that would be genuinely intuitive.  
The following answer repeats the original (and standard) development of the Poisson distribution as a limiting value of Binomial distributions, but instead of trying to obtain the probabilities directly (which requires some knowledge of the Gamma function), it focuses on a particularly simple relationship among the various probabilities.  The desired result drops out immediately with no calculation at all.

Imagine a large number $n$ of bins and a coin with a chance $p$ of falling heads.  Flip the coin independently for each bin, "filling" it when the coin falls heads and otherwise leaving it "empty."  Suppose $k$ of the bins end up full as shown with the dark shading in the figure, and the remaining $n-k$ are empty.

This is the situation depicted at the right (showing $k=4$ full bins).  The situation at the left is another possible outcome with only $k-1$ full bins.  It was created by choosing one of the $k$ bins at the right.  The "$k,1-p$" in the top arrow reminds us of the two chances involved in relating the right to the left: pick one bin out of $k$ (with chance $1/k$ for each such bin) and change the coin's outcome (with chance $1-p$).
The bottom arrow shows what changes in moving from the left to the right: one of the $n-(k-1)$ empty bins is chosen and then filled by an outcome having probability $p$.
This shows how to relate the chance of filling $k$ bins, which I will write $\pi(k,n,p)$, to the chance of filling $k-1$ bins, $\pi(k-1,n,p)$:
$$ \pi(k,n,p)(k)(1-p) = \pi(k-1,n,p)(n-k+1)(p).\tag{1}$$
Suppose now that $n$ becomes arbitrarily large, but as it does, $\lambda = p n$ (the expected proportion of full bins) stays constant.  Multiplying both sides by $n/k$ and dividing them by $n-k+1$ preserves the equality $(1)$, enabling us to rewrite it in terms of $\lambda$ with all the dependency on $n$ on one side:
$$\pi(k,n,n\lambda)\left(\frac{n-\lambda}{n-k+1}\right) = \pi(k-1,n,n\lambda)\left(\frac{\lambda}{k}\right).\tag{2}$$
For a fixed value of $k$, as $n$ grows the coefficient on the left of $(2)$ becomes uniformly close to $1$.  Thus, to an excellent approximation (which is on the order of $1/n$), we may take
$$\pi(k,n,n\lambda) \approx \pi(k-1,n,n\lambda)\left(\frac{\lambda}{k}\right)\tag{3}$$
for this particular $k$ (and all smaller values, too).
This, of course, recapitulates the construction of the Poisson distribution as a limit of Binomial$(n, \lambda/n)$ distributions.  We may now construct the full limiting distribution in terms of the chance of $k=0$ for the Poisson$(\lambda)$ distribution,
$$p_\lambda(0) = \lim_{n\to\infty}\pi(0,n,n\lambda).$$
According to $(3)$, this is multiplied by $\lambda/1$ to obtain $p_\lambda(1)$, and that is multiplied by $\lambda/2$ to obtain $p_\lambda(2)$, and so on.  Because this works for any finite $k$, it works for all $k$.

(The subscript $\lambda$ is dropped in the figure for brevity.)
For as long as $k \lt \lambda$, the chances keep increasing.  Once $\lambda \gt k$, they decrease.  Therefore the Poisson distribution is unimodal. 
In some special cases, the mode can occur at two adjacent values: this is precisely when $\lambda/k = 1$, for then the two chances $p_\lambda(k-1)$ and $p_\lambda(k)$ are equal.  Obviously this happens if and only if $\lambda$ is integral, in which case $k=\lambda$, QED.

Note that this is a stronger result than stated in the question, because it demonstrates the converse: when a Poisson distribution's mode occurs at two values, its rate $\lambda$ must be integral.
Incidentally, $p_\lambda(0)$--and therefore all the chances--is found by normalizing the sum 
$$\sum_{k=0}^\infty p_\lambda(k) = p_\lambda(0)\left[1 + \frac{\lambda}{1} + \left(\frac{\lambda}{1}\frac{\lambda}{2}\right) + \cdots + \left(\frac{\lambda}{1}\frac{\lambda}{2}\cdots\frac{\lambda}{k}\right)+\cdots\right]$$
to unity.  This sum is easily recognizable as $p_\lambda(0)\exp(\lambda)$, whence finally $p_\lambda(0) = 1/\exp(\lambda) = \exp(-\lambda)$ and therefore
$$p_\lambda(k) = \exp(-\lambda)\frac{\lambda^k}{1\cdot 2\cdots (k-1)\cdot k} = \frac{\lambda^k e^{-\lambda}}{k!}$$
for all $k \ge 0$.
A: Recall the density function for a Poisson RV is $e^{-\lambda}\lambda^x / x!$, you plugged $\lambda$ and $\lambda-1$ in for $x$. It is simple arithmetic to show it's true and that's enough intuition for me. The intuition is that Poisson RVs, for $\lambda \ge 1$, have convex density. If one allowed $x$ to take non-integral values, the curve would show a mode at $\lambda-0.5$. Proving this is easy, but it's simple to connect this with the idea that once $x$ is "forced" to take integral values, there are two modes interpolating this values.
curve(exp(4)*4^x/gamma(x+1), from=0, to=10)
abline(v=3.5)


