MinimalEquationTheta := proc(A::Matrix, v::posint, f:=y(x))
  local m, e, e2, B, i, j, k, d, q;
  m := LinearAlgebra:-RowDimension(A);
  q := -min(map(ldegree, A, x))-1;
  if q < 0 then q := 0 end if;
  e := Vector[row](m); e[v] := 1;
  B := Matrix(m);
  B[1] := e;
  for i from 2 to m do
    e := simplify(map(e2 -> x*diff(e2,x), e) + x*e.A);
    B[i] := e;
  end do;
  e := simplify(map(e2 -> x*diff(e2,x), e) + x*e.A);
  B := LinearAlgebra:-Transpose(B);
  d := LinearAlgebra:-Determinant(B);
  if d = 0 then error "non-cyclic" end if;
  return collect(x^(q*m*(m+1)/2) * (
    d*theta(f,x,m) - `+`(seq(
      LinearAlgebra:-Determinant(Matrix(m, (i,j) -> `if`(j=k, e[i], B[i,j])))
      * `if`(k > 1, theta(f,x,k-1), f),
      k = 1..m
    ))),
    [seq(theta(f,x,m-k), k=0..m-1), f, x]
  );
end proc:

_Trunc := proc(p,n)
  local d, r, pp;
  pp := collect(p, x);
  d := degree(pp, x);
  if d < n then
    return pp;
  end if;
  if pp::`+` then
    return select(t -> degree(t, x)<n, pp);
  else
    return 0;
  end if;
end proc:

MinimalTruncatedEquationTheta := proc(A::Matrix, v::posint, f:=y(x))
  local m, d, LT, C, T, r1, r2, r3, i, q;
  m := LinearAlgebra:-RowDimension(A);
  d := max(map(degree, A, x));
  LT := MinimalEquationTheta(A,v,f);
  C := [coeffs(LT, [y(x),seq(theta(y(x),x,i), i=1..m)], 'T')]; T := [T];
  LT := 0;
  q := -min(map(ldegree, A, x))-1;
  if q <= 0 then
    r1 := d + 2 - q*(m-2); # m
    r2 := r1;              # 2..m-1
    r3 := d + 2 - q*(m-1); # 0
  else
    r1 := d + 2 + q*(m+1);
    r2 := d + 2 + 2*q;
    r3 := d + 2 + q;
  end if;
  for i from 1 to nops(T) do
    if T[i] = f then
      LT := LT + (_Trunc(C[i],r3)+O(x^r3)) * y(x)
    elif op(-1,T[i]) = m then
      LT := LT + (_Trunc(C[i],r1)+O(x^r1)) * T[i]
    else
      LT := LT + (_Trunc(C[i],r2)+O(x^r2)) * T[i]
    end if;
  od;
  return LT;
end proc:

ReduceSolutions := proc(L::list)
  local R, i, s;
  R := copy(L);
  for i from nops(R) to 2 by -1 do
    s := eval(R[i-1], coeff(eval(R[i-1],O=0),x,ldegree(eval(R[i-1],O=0),x))=0);
    if s = R[i] then R[i] := 0 end if;
  od;
  return remove(e->e=0, R);
end proc:

LaurentPartialSolution := proc(A::Matrix, v::posint, d)
  local m, q, B, td, B1, p, alpha, k, LT, L, i, j, f, R;
  m := LinearAlgebra:-RowDimension(A);
  q := -min(map(ldegree, A, x)) - 1;
  B := copy(A);
  if nargs < 3 then
    td := max(map(degree, A))
  else
    td := d;
    B := map(e -> _Trunc(e,d), B);
  end if;
  if q <= 0 then
    R := td + 2 - q * (m - 2)
  else
    R := td + 2 + q
  end if;
  LT := MinimalEquationTheta(B, v);
  if min(map(ldegree, [coeffs(LT, [seq(theta(y(x),x,k),k=1..m),y(x)])], x))>=R then
    error "System is not well-defined w.r.t. unknown # %1", v
  end if;
  p := td;
  do
    p := p + 1;
    B1 := Matrix(m, (i,j) -> B[i,j] + `+`(seq(alpha[i,j,k]*x^k, k=td+1..p)));
    LT := MinimalTruncatedEquationTheta(B1, v);

    L := TruncatedSeries:-LaurentSolution(LT, y(x));
    if not L::list then return L end if;
    f := true;
    for i from 1 to nops(L) do
      j := ldegree(select(has, L[i], alpha), x);
      if j = infinity then
        f := false;
        break;
      end if;
      L[i] := remove(has, eval(L[i], O=0), alpha) + O(x^j);
    od;
    if f then break end if;
  od;
  return ReduceSolutions(L);
end proc: