Problem:
$P(i,-1) = 0$ for $i \ge 0$ and $P(i,0) = 1$ for $i \ge 0$
$P(0,j) = 0$ for $j \ge 1$
$P(i,j) = P(i-1,j-2) + c.P(i-1,j-1) + P(i-1,j)$ for all $i \ge 1$ and all $j \ge 1$
and $c$ is a constant.
One thing I must add, I am only interested in $P(i,i)$ . I am not interested in those cases where $i \ne j$ which is why writing it as
$P(n,n) = P(n-1,n-2) + c.P(n-1,n-1) + P(n-1,n)$ made more sense for me.
My attempts:
Since I am only interested in $P(i,i)$ My idea is to express as
$P(n,n) = g(c).P(n-1,n-1) + g'(c)$ so that we can convert the 2D recurrence into 1D.
Since first index i only depends of the earlier i-1 for all the three terms on RHS, so we should be able to get rid of it and convert it into 1D recurrence. Then it might be solvable using generating function.
Evaluating few elementary cases from bottom I find,
$P(1,1) = c$
$P(2,2) = 2 + c^2$
$P(3,3) = 2c+ c(2 + c^2)+ 2c = c(c^2+6) = 6c + c^3 $
$P(4,4) = (1+c(2c)+(2+c^2)) + c(6c + c^3) + ((2+c^2)+c(2c)+1) = 6 + 12c^2 + c^4$
$P(5,5) = P(4,3)+c.P(4,4)+P(4,5) = [P(3,1)+cP(3,2)+P(3,3)]+c.P(4,4)+[P(3,3)+cP(3,4)+P(3,5)] = P(3,1)+c.P(3,2)+2P(3,3)+c.P(3,4)+P(3,5) + c.P(4,4) = P(3,1) + c[P(2,0)+c.P(2,1)+P(2,2)] + 2P(3,3)+ c[P(2,2)+c.P(2,3)+P(2,4)] + (P(2,3)+cP(2,4)+P(2,5)) + c.P(4,4) = P(3,1) + c.P(2,0) + c^2.P(2,1) + c.P(2,2) + 2P(3,3) + c.P(2,2) + c^2.P(2,3) + c.P(2,4) + P(2,3) + cP(2,4) + P(2,5) + c.P(4,4) = P(3,1) + c.P(2,0) + c^2.P(2,1) + 2cP(2,2) + (c^2+1)P(2,3) + 2c.P(2,4) + P(2,5) + 2P(3,3)+cP(4,4) = 3c + c + c^2.3c + 2c.(2 + c^2) + (c^2+1).2c + 2c.1+ 0 + 2.(6c + c^3)+c(6 + 12c^2 + c^4) = 3c + c + 3c^3 + 4c + 2c^3 + 2c^3 + 2c + 2c + 12c + 2c^3 + 6c +12c^3 + c^5 = 30c + 21c^3 + c^5$
I tried to expand couple of steps manually from top as well, I am getting something like binomial coefficients as coefficients of the expanded terms.
Moreover we should be able to express it as
$P(n,n) = \sum_{k=1}^{n-1}g_k(c).P(n-k,n-k) + g'(c)$
I also tried programatically, but I don't know how to write P(i, 0):1 for all i>0 or how to write P(0,j):0 for all j>0 in the program.
from sympy import *
P = Function('P')
n = symbols('n',integer=True)
c = symbols('c', constant=True)
f = P(n,n) - P(n-1,n-2) - c*P(n-1,n-1) - P(n-1,n)
print(rsolve(f,P(k,k),{P(0,0):1, P(0,1):0, P(0,2):0}))
If possible please also do suggest a programatic way or fix my python code above.
Let $f(x,y)$ denote the [standard] generating function of $P(i,j)$. That is, $$f(x,y) = \sum_{i,j\geq 0} P(i,j)x^iy^j.$$
Using the definition of $P(i,j)$ to guide us, we manipulate this expression as follows: $$\begin{align*} f(x,y) &= \sum_{i\geq 0} P(i,0)x^i + \sum_{j\geq 1} P(0,j)y^j + \sum_{i,j\geq 1} P(i,j)x^iy^j \\ &= \sum_{i\geq 0}1\cdot x^i + \sum_{j\geq 1}0\cdot y^j + \sum_{i,j\geq 1}(P(i-1,j-2)+cP(i-1,j-1)+P(i-1,j))x^iy^j \\ &= \sum_{i\geq 0}x^i+xy^2\sum_{i\geq 0,j\geq -1}P(i,j)x^iy^j + cxy\sum_{i,j\geq 0}P(i,j)x^iy^j+x\sum_{i\geq 0,j\geq 1} P(i,j)x^iy^j \\ &= \sum_{i\geq 0}x^i+xy^2f(x,y)+cxyf(x,y) + x\left(f(x,y) - \sum_{i\geq 0}x^i\right) \\ &= 1 + (xy^2+cxy+x)f(x,y). \end{align*}$$
Solving for $f(x,y)$ yields $$f(x,y) = \frac{1}{1-xy(y+c+y^{-1})} = \sum_{i\geq 0}(xy(y+c+y^{-1}))^i$$ and therefore $$\begin{align*} P(n,n) &= [x^ny^n]\sum_{i\geq 0}x^iy^i(y+c+y^{-1})^i \\ &= [y^0](y+c+y^{-1})^n \\ &= [y^0]\sum_{k=0}^{\lfloor n/2\rfloor}\binom{n}{2k}(y+y^{-1})^{2k}c^{n-2k} \\ &= \sum_{k=0}^{\lfloor n/2\rfloor}\binom{n}{2k}\Big([y^0] (y+y^{-1})^{2k}\Big)c^{n-2k} \\ &= \sum_{k=0}^{\lfloor n/2\rfloor} \binom{n}{2k}\binom{2k}{k}c^{n-2k}. \end{align*}$$