Existence of Solution, System of Equations Suppose $P(\lambda, i)$ is the probability that a Poisson random variable with average $\lambda$ is equal to $i$, i.e. $\frac{\lambda^i}{e^{\lambda}i!}$
I think the following system of equations always has solution in $x$ and $y$, non-negative real numbers, for any $\alpha>0$ and $k\in \mathbb{N}_+$
\begin{cases}
\alpha=\sum_{i=0}^{\infty}P(x, i)\cdot P(y, k+i)  \\
\alpha=\sum_{i=0}^{\infty}P(x, i)\cdot P(y, k+i+1)
\end{cases}
where the necessary condition $\alpha\leq P(k+1, k+1)$ holds. It is easy to prove that this is indeed a necessary condition, equivalent to the condition that $\alpha=P(\lambda,k+1)$ has a solution. It is also easy to see that solution $y$ of the system is smaller or equal to $\lambda$, the largest solution to the equation $\alpha=P(\lambda,k+1)$. Experiments show that for each fixed $\alpha$ and $k$, there is a solution, but I did not manage to prove it analytically. 
Is there any analogue for mean value theorem for multidimensional functions?
Any suggestion for the proof directions will be appreciated.  
 A: Let $a:=\alpha$ and 
\begin{equation*}
 F_k(x,y):=\sum_{j=0}^\infty \frac{x^j}{j!}\frac{y^{k+j}}{(k+j)!}\,e^{-x-y},
\end{equation*}
assuming the standard convention $0^0:=1$. 
We have to consider the existence of a solution in $x$ and $y$ of the system 
\begin{equation*}
 a=F_k(x,y)=F_{k+1}(x,y). \tag{1}
\end{equation*}
We shall prove the following. 

Theorem 1. Take any natural $k$ and any 
  \begin{equation*}
 a\in(0,a_k],\quad\text{where}\quad a_k:=\sup_{x,y\ge0}F_{k+1}(x,y). \tag{1.5}
\end{equation*}
  Then the system (1) has a solution $x,y\ge0$. 

Remark 1. Since $F_k>0$, the condition $a\in(0,a_k]$ is obviously necessary in Theorem 1.
Proof of Theorem 1. Note that $F_k(x,y)\ge0$ for any real $x,y\ge0$ and $F_k(x,y)$ is continuous in real $x,y\ge0$. The crucial observation is the identity 
\begin{equation*}
 \partial_y F_{k+1}(x,y)=F_k(x,y)-F_{k+1}(x,y) \tag{2}
\end{equation*}
for real $x,y$. 
Next, fix for a moment any real $x\ge0$. Then $F_{k+1}(x,0)=0$ and, by dominated convergence, $F_{k+1}(x,\infty-)=0$. So, $F_{k+1}(x,y)$ attains its maximum in $y$ at some real point $y=y_x\ge0$. At this point, we have $\partial_y F_{k+1}(x,y)=0$. So, by (2), 
\begin{equation*}
 F_k(x,y_x)=F_{k+1}(x,y_x)=\max_{y\ge0}F_{k+1}(x,y)=:M_{k+1}(x), \tag{3} 
\end{equation*}
for all real $x\ge0$. 
Next, 
\begin{align*}
 M_k(x)&\le\sum_{j=0}^\infty \frac{x^j}{j!}\max_{y\ge0}\frac{y^{k+j}}{(k+j)!}\,e^{-x-y} \\ 
 &=\sum_{j=0}^\infty \frac{x^j}{j!}e^{-x}\frac{(k+j)^{k+j}}{(k+j)!}\,e^{-k-j} \tag{4} \\ 
 &\ll\sum_{j=0}^\infty \frac{x^j}{j!}e^{-x}\frac1{\sqrt{k+j}}=E\frac1{\sqrt{k+\Pi_x}}
 \underset{x\to\infty}\longrightarrow0 
\end{align*}
by dominated convergence and because $\Pi_x\underset{x\to\infty}\longrightarrow\infty$ in probability, where $\Pi_x$ is a Poisson random variable with parameter $x$. So,
$M_k(\infty-)=0$. It is also not hard to see that 
$F_k(x,y)$ is continuous in real $x\ge0$ uniformly in real $y\ge0$ (see the Appendix), so that $M_k(x)$ is continuous in $x\ge0$. So, $M_{k+1}(x)$ attains its maximum in $x\ge0$ (equal $a_k$, by (1.5)) and takes all values in the interval $(0,a_k]$. Now Theorem 1 follows by (3). $\qquad\Box$
Appendix. Similarly to (2), 
\begin{equation*}
 \partial_x F_k(x,y)=F_{k+1}(x,y)-F_k(x,y). 
\end{equation*}
for real $x,y$. Therefore and because $0\le F_k\le1$, we have $|\partial_x F_k(x,y)|\le1$ for real $x,y$, so that $F_k(x,y)$ is indeed continuous in real $x\ge0$ uniformly in real $y\ge0$. 
Added: Let us now show that 
\begin{equation*}
 a_k=c_{k+1},\quad\text{where}\quad c_k:=\frac{k^k}{k!}\,e^{-k}
 \sim\frac1{\sqrt{2\pi k}}   
\end{equation*}
as $k\to\infty$. To this end, note first that 
\begin{equation*}
 c_{k+1}/c_k=(1+1/k)^k/e<1,
\end{equation*}
and so, $c_k$ is decreasing in $k$. 
So, recalling (4), we have 
\begin{equation*}
 M_k(x)\le\sum_{j=0}^\infty \frac{x^j}{j!}e^{-x}\,c_{k+j}
 \le\sum_{j=0}^\infty \frac{x^j}{j!}e^{-x}\,c_k=c_k=M_k(0). 
\end{equation*}
Thus, in view of (1.5) and (2), 
\begin{equation*}
 a_k=\max_{x\ge0}M_{k+1}(x)=M_{k+1}(0)=c_{k+1}, 
\end{equation*}
as desired. 
In particular, for $k=0,1,2,3$ the values of $a_k$ are $\approx0.367879, 0.270671, 0.224042$. 
A: The sums over $P(\lambda,i)=\frac{\lambda^i}{e^{\lambda}i!}$ are evaluated in terms of a Bessel function as
$$\sum_{i=0}^{\infty}P(x, i)\cdot P(y, k+i)=y^k e^{-x-y} \left(xy\right)^{-k/2} I_k\left(2 \sqrt{xy}\right)$$
$$\sum_{i=0}^{\infty}P(x, i)\cdot P(y, k+i+1)=(y/x)^{1/2}\,y^{k} e^{-x-y} \left(xy\right)^{-k/2} I_{k+1}\left(2 \sqrt{xy}\right)$$
for any positive integer $k$ these two expressions should be equal for some $x,y>0$. (For $x=y=0$ both expressions are identically zero.) So the function
$$F_k(x,y)=\sqrt{x}\, I_k\left(2 \sqrt{xy}\right)-\sqrt{y} \,I_{k+1}\left(2 \sqrt{xy}\right)$$
should pass through zero in the quadrant $x,y>0$ for any positive integer $k$.
For large $z=2\sqrt{xy}$ both Bessel functions $I_k(z)$ and $I_{k+1}(z)$ grow as $(2\pi z)^{-1/2}e^z$, so by making $x$ much larger than $y$ the function $F_k(x,y)$ is positive and by making $y$ much larger than $x$ it is negative, hence it must go through zero when $x\approx y$.

I had not appreciated that $\alpha$ is fixed from the beginning like $k$, not a variable like $x$ and $y$. So we also need to show that $x\approx y\gg 1$ allows the sum to reach any $\alpha>0$, so
$$\alpha=e^{-2x}  I_k\left(2x\right)\approx (4\pi x)^{-1/2},\;\;x\gg 1.$$
This is possible only for $\alpha\ll 1$. The OP lists as necessary condition $\alpha\leq P(k+1,k+1)$, it is not clear to me this is sufficient.
