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.

Could you please help me? Thank you very much in advance!

A minimal example of A(z), psi(z) and p could be: (x is an external constant)

n=3;
A(z_) = {{z, 1/z, 1/(z - 1)}, {z^2, x z, (1 - x) z}, {0, 1/(x z), z}};
psi(z_) = {1, 1/z, 1/(z-1)};
p=1;