Let us suppose that I have some differential operator $D_{x,y}$ acting on functions of two variables $(x,y)$. I want to solve one eigenvalue equation
$$D_{x,y}f_{ij}(x,y)=\lambda_{i,j} f_{ij}(x,y).$$
I know that the $\lambda_{ij}$ and I know that $f_{ij}$ admit a power series expansion
$$f_{ij}(x,y)=\sum_{m,n=0}^\infty a_{mn}x^{\alpha+m}y^{\beta+n}.$$
I need to find the coefficients $a_{mn}$, which everyone in the literature says its very easy to do using the differential equation with a computer. I want to know how to do that using Mathematica, because it feels to me that it should be possible. Even by hand it is possible, it is just that it is very tedious and I wanted to use Mathematica for this.
What is the way to use Mathematica in order to find the coefficients of a power series expansion of a function that solves some differential equation? How can I find the $a_{mn}$, $\alpha$ and $\beta$, given $D_{x,y}$ and $\lambda_{ij}$?
AsymptoticDSolveValue, but I don't think it works on PDEs. Another possibility: If you equate and truncate two series, you can solve for the coefficients withSolveAlways. Or equate theCoefficientList[]of each. – Michael E2 Jun 23 '23 at 17:56AsymptoticDSolveValueapproach for simple ODEs: https://mathematica.stackexchange.com/a/96696/4999. I avoidedSolveAlwaysbecauseSolveAlwayscan choke on a large system – Michael E2 Jun 23 '23 at 18:02