Is this idea for solving multivariable recurrences close to correct? According to Wikipedia, the binomial coefficient identity offers a good example of a multivariable recurrence relation.  I believe the idea is to treat the top and bottom slots of the choose operator each as an input to an unknown function, leading to the equation $$C[n,\ k] = C[n - 1,\ k - 1] + C[n - 1,\ k]$$
where we pretend not to know in advance that $C$ is the choose operator.  My goal is to try to solve this analytically by extrapolating from the patterns relating discrete and continuous equations in the single variable case.
First, I shifted some things to get $C[n + 1,\ k + 1] - C[n,\ k + 1] - C[n,\ k] = 0$, and considered the most similar partial differential equation, $c_{_{NK}} - c_{_K} - c = 0$, where $c$, $N$, and $K$ are the continuous versions of $C$, $n$, and $k$, respectively.  This was based on the observation that incrementing inputs and taking derivatives have the same effect on the characteristic equation in linear, single variable equations.
Next, I tried to solve the partial differential equation using separation of variables.  This led to $c = (Ae^{\frac{1}{\lambda} N})(Be^{\lambda K + K})$, where $A$ and $B$ are the arbitrary constants from each separated ODE.
Finally, I want to get rid of the exponential nature of the solution and write the solution to the discrete equation.  The most natural thing for me to write based on the single variable pattern is $$C = a(\frac{1}{\Lambda})^n b(\Lambda + 1)^k$$
I just yeeted the $e$ from each factor and pulled the coefficients of the independent variables down (and switched the case of each letter again for clarity of notation).  I'm not sure if $a$ and $b$ can be lumped, since they may have been affected by two different discretizations.
Actually, I'm not sure if any of this is valid.  Is this process or something similar a valid way to solve multivariable recurrence relations?  If not, why is there a nice relationship between analytically solving discrete and continuous equations in the single variable case but not in the multivariable case?
 A: You can use the wonderful theory of generating functions to make this all rigorous! Let
$$
c(x,y)=\sum_{n=0}^\infty \sum_{k=0}^\infty C[n,k]\frac{x^n}{n!}\frac{y^k}{k!}
$$
This is standard practice; $c(x,y)$ is known as the bivariate exponential generating function of the bivariate sequence $C[n,k]$. I will write this as
$$
c(x,y)\newcommand{\gen}{\quad \text{generates}\quad }\gen C[n,k]
$$
You can then show, by manipulating these power series, that
$$
\begin{align}
c(x,y)\gen C[n,k]
&\iff
c_x(x,y)\gen C[n+1,k]
\\&\iff
c_y(x,y)\gen C[n,k+1] 
\end{align}
$$
Now, consider the following function:
$$
\sum_{n=0}^\infty \sum_{k=0}^\infty \big(C[n+1,k+1]-C[n,k+1]-C[n,k]\big)x^ny^k
$$
On the one had, the recurrence relation you have implies that this function is zero. On the other hand, we can split this into three sums, each of which is a certain partial derivate of $c$. The result is
$$
c_{xy}-c_y-c=0,\tag{*}
$$
exactly the equation you arrived at. As you did, we can try assuming that $c(x,y)=a(x)b(y)$, in which case this becomes
$$
a'b'-ab'-ab=0,\\
a'b'-ab'+\color{blue}{a'b}-ab-\color{blue}{a'b}=0,\\
(a'-a)(b'+b)-a'b=0,\\
\frac{a'-a}{a'}=\lambda=\frac{b}{b'+b},\\
a'=\frac{1}{1-\lambda}a,\qquad b'=(\tfrac{1-\lambda}{\lambda}) b,\\
a(x)=\exp(\tfrac{1}{1-\lambda}x)\qquad b(y)=\exp(\tfrac{1-\lambda}{\lambda}y)
$$
Now that we have solved for $c(x,y)$, we need to recover $C[n,k]$. This is done by extracting the coefficient of $x^n$ in $a(x)$, and multiplying that with coefficeint fo $y^k$ in $b(y)$. The result is
$$
C[n,k]\stackrel{?}= (1-\lambda)^{k-n}\lambda^{-k}
$$
For example, if you choose $\lambda=\frac12$, then you get the solution $C[n,k]=2^n$. This certainly satisfies the equation $C[n+1,k+1]-C[n,k+1]-C[n,k]$.
However no value of lambda gives the correct solution $C[n,k]=n!/(k!(n-k)!)$, so what went wrong? The problem was when we assumed that the solution to $(*)$ was separable. In fact, assuming that $c(x,y)=a(x)b(y)$ will only generate solutions to $C[n,k]$ of the form $A[n]B[k]$. Since $n!/(k!(n-k)!)$ is not of this form, you will need fancier PDE methods to actually be able to solve this equation. Specifically, you will need to use the boundary conditions
$$
c(x,0)
=\sum_{n=0}^\infty C[n,0]\frac{x^n}{n!}
=\sum_{n=0}^\infty 1\cdot \frac{x^n}{n!}=e^x,\\
c(0,y)=\sum_{k=0}^\infty C[0,k]y^{k}{k!}=1+0y+0y^2/2!+\dots=1$$
in order to determine the correct solution to $(*)$. Unfortunately, I do not know how to proceed here. But this general method is certainly useful in many situations, and provides a good bridge between the continuous and the discrete.
