Peak of Poisson distribution I was experimenting with Poisson distribution on desmos and noticed that as $\lambda$ (the mean of distribution) gets larger, the x coordinate of distribution peak approaches $\lambda-0.5$. Why is it so?

The analysis at https://stats.stackexchange.com/a/211612/919 shows the mode of any Poisson$(\lambda)$ distribution is near $\lambda$ itself.  Although that question concerns only integral $\lambda,$ its results answer the present question, too.
Let $p_\lambda(k)$ be the Poisson$(\lambda)$ probability for $k\in\{0,1,2,\ldots,\}$ (all other probabilities are zero, of course).  These probabilities are proportional to the ratios
$$p_\lambda(k) \ \propto\  \frac{\lambda^k}{k!}.$$
As pointed out in the foregoing link, this means two successive Poisson probabilities are related by
$$p_\lambda(k+1) = \frac{\lambda}{k+1} p_\lambda(k).$$
Consequently, as $k$ progresses from $0$ through all integers less than $\lambda-1,$ the probability increases from $p_\lambda(k)$ to $p_\lambda(k+1);$ and then once $k$ exceeds $\lambda-1,$ the probability decreases.  Therefore

All Poisson probability functions $p_\lambda$ rise to a peak and then fall again.  The peak (a mode) occurs at the greatest integer less than $\lambda,$ written $\lfloor \lambda\rfloor.$  When $\lambda$ is an integer, the peak occurs at the two neighboring values $\lfloor \lambda\rfloor$ and $\lfloor \lambda\rfloor + 1.$

To illustrate, here is a plot of the modes based on a brute-force search.

The search was conducted by this R function.  It evaluates the Poisson probabilities for $k$ between two suitable extreme quantiles, finds the indexes where the largest probability occurs, and returns an array of the values of $k$ corresponding to those indexes.
mode <- function(lambda) {
  q <- min(1/2, ppois(floor(lambda), lambda))          # This is *a* probability
  k <- do.call(seq, as.list(qpois(c(q, 1-q), lambda))) # Search limits
  p <- dpois(k, lambda)
  k[which(p == max(p))]
}

