Algorithm for using power series to numerically solve a partial differential equation given a boundary condition?
Instead of building a polynomial and trying to deduce the equation on the coefficients, you can just derive it directly from a general term coefficient

&gt;&gt;&gt; i, j, x1, x2 = symbols('i j x1 x2')
&gt;&gt;&gt; a = IndexedBase('a')
&gt;&gt;&gt; u = Function('u')
&gt;&gt;&gt; pde = u(x1, x2).diff(x1, 2) - u(x1, x2).diff(x2, 2)
&gt;&gt;&gt; term = a[i, j]*x1**i*x2**j
&gt;&gt;&gt; pdeterm = powsimp(pde.replace(u, Lambda((x1, x2), term)).doit())
&gt;&gt;&gt; pdeterm
i*x1**(i - 2)*x2**j*(i - 1)*a[i, j] - j*x1**i*x2**(j - 2)*(j - 1)*a[i, j]


The PDE being identically zero means each coefficient of each x1, x2 term is 0. 

&gt;&gt;&gt; pdeterm.coeff(x1**i*x2**(j - 2))
-j*(j - 1)*a[i, j]
&gt;&gt;&gt; pdeterm.coeff(x1**(i - 2)*x2**j)
i*(i - 1)*a[i, j]


This is the same as what you had above if you shift i and j

&gt;&gt;&gt; pdeterm.coeff(x1**i*x2**(j - 2)).subs({j: j + 1})
-j*(j + 1)*a[i, j + 1]
&gt;&gt;&gt; pdeterm.coeff(x1**(i - 2)*x2**j).subs({i: i + 2})
(i + 1)*(i + 2)*a[i + 2, j]


I'm not sure if there's an advantage to having them one way or another.  For a general PDE, you'd want to write a function that does this automatically.

The boundary conditions can be managed in a similar way using subs. 

You can then solve the system analytically with solve or numerically with nsolve.