# differential equations – Local Series solution at singular point for system of first order ODEs

I want to calculate `Psi(z)` in the equation

``````D(Psi(z), z) + A(z).Psi(z) == psi(z)
``````

around a given point `p`, where `A(z)` is an $$ntimes n$$–matrix, `psi(z)` is an $$n$$–vector and hence `Psi(z)` also needs to be an $$n$$–vector. The matrix `A(z)` and vector `phi(z)` can both be singular at `p`, but neither has an essential singularity. Can you please help me solve this problem efficiently?

The goal is then to take `Psi(z)` and multiply it with some given `phi(z)`, and then take the residue at `p`. Therefore ideally the outcome is a `SeriesData` object around `p`, or something similar to that. Note that `p` generally is a point on $$CP^1$$, so `p` can be `ComplexInfinity`.

I have tried a few things and some were succesfull in subcases, but while those often worked in reasonable time for $$n=1$$ they took unreasonably long for $$n>4$$ or had similar issues.

• I generally wrote `A(z)` as a list of `SeriesData` objects with coefficients `Psico(l, k)` $$(l in {1, dots, n})$$ and expanded `A(z)` and `psi(z)` using `Series` and took the list of their coefficients `Acoef` and `psiCoef`, which reduces the equation to something of the form
``````(a+1) Psico(j, a+1) + Sum(Acoef(b)((j, k)) * Psico(k, a-b), {b, low, high}, {k, 1, n}) == psiCoef(j, a).
``````
• For `A(z)` non-singular in `z`, the equation becomes an easy recursive equation in the expansion coefficients of `Psi(z)`. This can be solved either using recursive methods (such as `RecurrenceTable` or a `For` loop by hand) or by writing the equations as a triangular matrix and solving that (using e.g. `LinearSolve`). However, this ran into problems when `A(z)` is singular at `p` as in that case the matrix is no longer triangular and equivalently the recursive method requires coefficients which we do not yet know.
• For $$n$$ small, the matrix could still forcefully be solved using `LinearSolve`, but for large $$n$$ that becomes incredibly slow. Perhaps that is because my code was too naive, and there are some tricks to convince LinearSolve that the matrix is nearly triangular?
• I then tried using methods such as `AsymptoticDSolveValue`, but that persistently tried to approximate `Psi(z)` with a polynomial (not using singular terms). I am using Mathematica 12.0 Student Edition (where `AsymptoticDSolveValue` is `Experimental`); I have understood this method has been changed slightly in the newer version so perhaps this is now better but I have no acces to newer version of Mathematica.

Finally remark that there is of course not just a single solution, but your can generally add the solution to the homogeneous equation. This can be ignored, as it is mathematically expected that the homogeneous solution does not contribute to the later residue. Therefore not all solutions are needed, but just one solution should be sufficient. Of course, the homogeneous solution generally has an essential singularity, so that solution can typically be avoided somewhat by choosing some sufficiently low coefficient `Psico(k, low)` and fixing it to be 0 by hand.

A minimal example of `A(z)`, `psi(z)` and `p` could be: (`x` is an external constant)
``````n=3;