Algorithm for using power series to numerically solve a partial differential equation given a boundary condition? Motivation:
Following this discussion about using asymptotic expansions (i.e. polynomial power series) for numerically solving partial differential and algebraic equations (PDAE), I couldn't find any implementation of the method. So I'm thinking of implementing a SymPy function similar to Mathematica's AsymptoticDSolveValue (here)  but for PDEs (here or here). So far I'm able to generate a symbolic multivariate polynomial given a list of non-negative integers D=[d1,...,dm] (here, here and here).
Example:
Now I'm able to use the symbolic multivariate polynomial to numerically solve a PDE. For example (from here) given the PDE:
$$\frac{\partial^2 u}{\partial x_1^2} - \frac{\partial^2 u}{\partial x_2^2}=0 \, ,\tag{1}$$
and the boundary conditions:

$u(x_1,0)=x_1^2+x_1\, , \tag{2}$
$u_{x_2}(x_1,0)=2x_1+1 \, , \tag{3}$



I can generate a 2D symbolic polynomial in Sympy:
from itertools import product
from sympy import IndexedBase, symbols, Poly
D = (5, 5) # 5 and 5 are just some random integers, any non-negative integer should do
d = len(D)
indices=list(product(*map(range, D)))
a = IndexedBase('a')
coeffs = {i: a[i] for i in indices}
vars = symbols(f'x1:{d+1}')
u = Poly(coeffs, *vars)

Equation 1 can be implemented as:
pde=u.diff(0,0).add(-u.diff(1,1))

Implying
$$j(j+1)a_{i,j+1}=(i+1)(i+2)a_{i+2,j} \tag{4}$$
The first boundary condition:
u.eval(1,0)

Implying
$$a_{0,0}=0 \, , \, a_{1,0}=1 \, , \, a_{2,0}=1 \, , \,a_{i,0}=0 \, \forall \,(2&lt;i) \tag{5}$$
and the second boundary condition
u.diff(1).eval(1,0)

giving
$$a_{0,1}=1 \, , \, a_{1,1}=2 \, , \, a_{i,1}=0 \, \forall \, (2&lt;i) \tag{6}$$
Now from this point it is just a system of nonlinear algebraic equation of $a$s (in this specific case linear). Which should be solvable with other analytical/numerical methods.
Question:
I want to automate the process above. I want to have a function:
AsymptoticPdeSolve(eqns,fs,vars,D)

Where eqns is the set of symbolic partial differential expressions, fs are the set of functions we want to solve, vars are the variables and D is the dimension of our multivariate polynomials.
I would appreciate if you could help me know what is the best algorithm for this process. I will use an existing symbolic library/software like SymPy, so anything already existing doesn't need to be reimplemented.
P.S. Axiom-FriCAS has a seriesSolve function which source code can be found here. Also since I posted this question, Nicolas CELLIER was so kind to implement an early version in sympy which can be seen in this Jupyter notebook.
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.