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. 



Existence of Solution, System of Equations
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$. 


