The hard part of the question is the second part about exhibiting a non-stationary, weakly stationary series with identical distributions at all times.
I hope you won't mind a simplification: rather than using Poisson distributions, which can take on infinitely many values, let's use one that takes on only a small number of values. A particularly simple distribution $F$ describes the difference between the number of heads and number of tails in two independent throws of a fair coin. (This makes it a version of a Binomial distribution.) It therefore has a chance of $1/2$ of being $0$ and a $1/4$ chance each of being $-2$ and $2$. You can readily compute that its expectation is $0$ and its variance is $2$.
I need to describe two different bivariate distributions, $G$ and $H$, whose marginals both equal $F$.
$G$ is the joint distribution of two independent random variables $X,Y$ distributed according to $F$. It assigns positive probabilities to all nine distinct outcomes $(X,Y)$: when both of $|X|$ and $|Y|$ are $2$, the chance is $1/16$; the chance of $(0,0)$ is $1/4$; and the other four chances are $1/8$.
$H$ assigns positive probabilities only to five of these outcomes: when both of $|X|$ and $|Y|$ are $2$, the chance is $1/8$, and the chance of $(0,0)$ is $1/2$.
(One way to realize $H$ in a simulation is to draw $X$ from $F$, and then negate $X$ with probability $1/2$ to produce $Y$.)
It is immediate that the covariance of $X$ and $Y$ under the law $G$ is zero, because $X$ and $Y$ are independent. Under the law $H$, the covariance also is zero, because (since the expectations are zero), $$\operatorname{Cov}_H(X,Y)=\mathbb{E}_H(XY) - \mathbb{E}_H(X)\mathbb{E}_H(Y) = \mathbb{E}_H(XY)$$
and a direct calculation from the definition of expectation gives $$\mathbb{E}_H(XY) = \frac{1}{8}(2\times 2 + 2\times (-2) + (-2)\times 2 + (-2)\times (-2)) + \frac{1}{2}(0\times 0) = 0.$$
We have established that $G$ and $H$ have the same first and second bivariate moments and identical marginal distributions. They provide a nice example of a weakly stationary time series $(Z_t),\ t\in\mathbb{Z},$ that is not stationary. Simply let $(X_t),\ t\in\mathbb{Z}\setminus\{0,1\}$ be a sequence of independent random variables with distribution $F$ and let $(X_0,X_1)$ be independent of the other $X_t$ and have joint distribution $H$. Thus, all bivariate distributions of $(X_t, X_s)$ follow the law $G$ except when $\{t,s\}=\{0,1\}$, which follow the law $H$. Since $H$ differs from $G$, this makes the time series non-stationary. (The illustrations at the end help you see this.)
Nevertheless, all the univariate distributions are the same and all the covariances between $X_t$ and $X_s$ for $t\ne s$ are zero, as we have seen. That makes the series weakly stationary.
One way to illustrate the situation is to show some partial realizations of this process. What is particularly revealing are the first differences of the time series, $dX_t = X_{t+1} - X_t$. These will be values in the set $\{-4,-2,0,2,4\}$. The construction of $H$ precludes the values $\pm 2$ from ever occurring for $dX_0 = X_1 - X_0$, but those values can occur for all other differences. In the graphics below I have therefore plotted these first differences against time $t=-15, -14, \ldots, 15$. The value $dX_0$ is highlighted in red. The "odd" values $\pm 2$ are plotted with crosses and the other values with dots.
It helps to see many realizations at once. In the next graphic I overlaid 300 of them on the same axes.
The pattern becomes clear: at time $t=0$, there are never any differences of $2$ or $-2$ plotted, but such differences show up at all other times. That's a visible demonstration of non-stationarity of the differences. Since the differences are not stationary, the original series cannot be either.
The R code to produce these simulations is particularly simple, so it might be useful to see it.
times <- (-15):16
n <- length(times)
n.reps <- 300
realizations <- replicate(n.reps, {
x <- (rbinom(n, 2, 1/2)-1)*2 # Draw independent Binomial variates
names(x) <- times
x["1"] <- sample(c(-1,1), 1)*x["0"] # Create a dependence at times 0 and 1
x
})
