FormalSolution := proc()
local check_arguments, EG, Row_Length, Shift_Row, ResolvingSequence, ResolvingStep, 
      input_type, ModuleApply;
# Input list/set form:
#    Sys  -- homogeneous ordinary differential system
#    vars -- list of unknown functions, [y1(x),y2(x)] for example
#   'select_indicator' = si -- (option) si is
#                              a name from vars or
#                              a procedure to select an unknown from vars,
#                              with arguments (system-in-set-form, list of vars)
# Output: a basis of formal solutions of the system Sys
# Input (matrix-form):
#    Sys  -- linear homogeneous equation with matrix coefficients
#       A_n . diff(y(x),x$n) + ... A_1 . diff(y(x),x) + A_0 . y(x) = 0
#    y(x) -- unknown function
#    'select_indicator' = si -- (option) si is
#                             an integer or
#                             a procedure to select integer,
#                             with arguments (system-in-matrix-form, var)
# Input (normal form):
#    A  -- square matrix
#    x  -- independent variable
#    'select_indicator' = si -- (option) si is an integer,
#                                   or a procedure to select integer by
#                                      the given matrix and x
#    'by_cyclic_vector' = true/false -- (option) apply a cyclic vector method,
#                                       false by default
# Output: a basis of formal solutions of the system y'(x) = A y(x)
#  where y is a vector of unknown functions
ModuleApply := proc()
local pp, st_total, st, L, ExpParts, ExpPart, t, t1, Res, new_vars, i, j,
      system, sys2, m, a, r, Q, lambda, ell, indts, slv, RegPart,
      constructed, si, terms, x, j1;
global _c;

    st_total := time();
    pp, t1, terms, i, si := check_arguments(args);
    x := pp['var'];

    st := time():
    L := ResolvingSequence(pp, i, si);
    userinfo(3, FormalSolution, "ResolvingSequence is constructed:", time()-st);

    Res := NULL;
    m := nops(pp['functions']);
    new_vars := [seq(z[j](t), j = 1 .. m)];
    system := pp['explicit_form'][0].Vector(pp['functions']) +
          seq(pp['explicit_form'][j].Vector(diff(pp['functions'],x$j)),
                                                          j = 1 .. pp['order']);
    system := convert(system, 'set');
    constructed := NULL;
    for ell in L do
        st := time():
        ExpParts := DEtools['gen_exp'](ell,
                  op(map2(op, 0, indets(ell, 'function')) intersect
                  {map2(op, 0, pp['functions'])[]})(x), t);
        userinfo(3, FormalSolution, nops(ExpParts),
                         " exponential parts is found:",
                                                     time()-st);


        ExpParts := map(el -> [lcoeff(lhs(el[-1]), t), degree(lhs(el[-1]), t),
                       el[1]-coeff(el[1], t, 0), coeff(el[1], t, 0)], ExpParts);
        ExpParts := map(el ->
                  `if`((el[4])::'rational', [el[1 .. 3][], frac(el[4])], el),
                                                                      ExpParts);
        for ExpPart in ExpParts do
            a, r, Q, lambda := ExpPart[];
            if [a, r, Q] in [constructed] then
                next;
            end if;
            constructed := constructed, [a, r, Q];
            ExpPart := int(r*Q/t, t);
            ExpPart := map(radnormal, ExpPart);
            userinfo(3, FormalSolution, "exponential part: ",
                                               [exp(ExpPart), x=a*t^r]);
            st := time():
            sys2 := PDEtools:-dchange({x = a*t^r,
                 seq(pp['functions'][j] = t^(r*lambda)*exp(ExpPart)*new_vars[j],
                                                      j = 1 .. m)}, {system[]});
            sys2 := radnormal(map(`*`, sys2, 1/(t^(r*lambda)*exp(ExpPart))));
            RegPart := LinearFunctionalSystems:-RegularSolution(
                                                            [sys2[]], new_vars);

            if {RegPart[]} = {0} then
                userinfo(3, FormalSolution, "no regular solution:", time()-st);
                next;
            end if;

            if terms > 0 then
                RegPart := LinearFunctionalSystems:-ExtendRegularSolution(
                                                              RegPart, terms-1);
            end if;
            if t1 = FAIL then
                slv := {seq(RootOf(a*t^r = x, t, 'index' = j),
                                                                     j = 1..r)};
                a := [seq(eval(t^(r*lambda)*exp(ExpPart), t = j).
                          eval(Vector(RegPart), t=j), j = slv)];
                indts := indets(a, 'O'('anything'));
                a := eval(a, map(el -> el = series(el, x), indts));
                indts := indets(a, series);
                a := eval(a, map(el -> el = O(x^op(2, el)), indts));
            else
                a := [t^(r*lambda)*exp(ExpPart).Vector(RegPart),
                      a*t^r = x];
            end if;

            indts := indets(a, _c['integer']);
            a := [seq(eval(eval(a, {j = 1}), {seq(j1 = 0,j1 = indts)})[],
                                                                    j = indts)];
            a := convert(a, 'radical');
            userinfo(3, FormalSolution, nops(a),
              " dimension regular solution space is found:", time()-st);

            Res := Res, a[];
        end do;
    end do;

    userinfo(3, FormalSolution, "all time:", time()-st_total);

    return [Res];
end proc:
ResolvingSequence := proc(pp, i, si)
local eq, B, m, n, d, st, vars, k, j, sys, slv, indts, S, A, tmp;

    m := nops(pp['functions']);
    n := pp['order'];

    eq, B := ResolvingStep(pp, i, m, n);

    d := LinearAlgebra:-RowDimension(B);
    if d = n*m+1 then
        return [eq=0];
    end if;

    # reduction step
    # solve B(x) Y(x) = 0 with yi(x)=0, ξyi(x)=0, ξ^{n-1}yi(x) = 0
    st := time():

    B := LinearAlgebra:-DeleteColumn(B[n+1..-2], [seq(i+m*k,k=0..n-1)]);
    vars := map2(op, 0, subsop(i = NULL, pp['functions']));
    vars := [seq(map(el-> el[k], vars)[], k = 0 .. n-1)];

    sys := convert(B.Vector(vars), 'set');
    slv := solve(sys, {vars[]});
    indts := indets(map(rhs, slv));
    indts := select(has, vars, indts);

    slv := remove(evalb, slv);
    slv := [seq(eval(vars[k], slv), k = 1 .. nops(vars))];
    S := Matrix([seq([seq(coeff(slv[k], indts[j]), j = 1 .. nops(indts))],
        k = 1 .. nops(vars))]);

    A := Matrix(convert(pp['normal_form'], 'list'));
    A := LinearAlgebra:-DeleteColumn(
           LinearAlgebra:-DeleteRow(A, i), [seq(i+m*k,k=0..n-1)]);
    A := LinearAlgebra:-SubMatrix(A,
             map(proc(el) local k;
               `if`(member(el, vars[-m+1..-1], 'k'), k, NULL) end proc, indts),
             [1..-1]).S;
    A := map(radnormal, A);

    A := <<LinearAlgebra:-SubMatrix(S,
       map(proc(el) local k;
          `if`(member(el, vars[1..-m], 'k'), k+m-1, NULL) end proc, indts),
           [1..-1]), A>>;

    tmp := check_arguments(A, pp['var']);
    userinfo(3, FormalSolution, "reduction is evaluated:", time()-st);
    tmp[1]['functions'] := map(el-> el(pp['var']), indts);

    if si::'integer' or si::'name' then
        eq := eq = 0, procname(tmp[1], 1, 1)[];
    elif input_type = 'Matrix' then
        si(A, pp['ring']);
        error "no tests";
    else
        error "no tests";
    end if;

    indts := indets([eq[2..-1]], 'function');
    indts := select(el-> op(0, el) in vars, indts);
    eq := eval([eq], map(el->op(0,el)=op([0,0],el), indts));
    if select(has, indets([eq[2..-1]]), vars[m..-1]) <> {} then
        error "no tests";
    end if;

    return eq;
end proc:
ResolvingStep := proc(pp, i, m, n)
local AA, k, B, slv, eq, st, C, indts;
eval(pp);
    userinfo(3, FormalSolution, "resolving for ", [0$(i-1), 1, 0$(m-i)]);

    AA := pp['normal_form'];

    st := time():

    B := Matrix(n,m*n,{seq((k+1,m*k+i)=1,k=0..n-1)});

    B := <<B, Matrix(1, m*n, [seq(AA[k][i], k=0..n-1)])>>;

    slv := LinearAlgebra:-LinearSolve(B^%T, Vector(m*n), 'free' = C);

    while convert(slv, 'set') = {0} do

       B := <<B,
              map(expand@radnormal, Matrix(1, m*n,
                 [B[-1,-m .. -1].AA[0] +
                  map(diff, B[-1,1 .. m], pp['var']),
        seq(
          B[-1,(k-1)*m + 1 .. k*m] +
          B[-1,-m .. -1].AA[k]+
          map(diff, B[-1,k*m + 1 .. (k+1)*m], pp['var']),
        k = 1..n-1)])) >>;
        slv := LinearAlgebra:-LinearSolve(B^%T, Vector(m*n), 'free' = C);

    end do;

    userinfo(3, FormalSolution, LinearAlgebra:-RowDimension(slv)-1,
                              "order resolving equation is constructed:",
                                                                     time()-st);

    indts := [indets(slv, C['integer'])[], C[0]];
    slv := eval(slv, {indts[1] = 1, seq(k = 0, k = indts[2..-1])});
    slv := map(expand@radnormal,
                    lcm(convert(map(denom, map(radnormal, slv)), 'set')[])*slv);

    eq :=  OreTools:-Converters:-FromOrePolyToLinearEquation(
               OrePoly(convert(slv, 'list')[]),
               op(0,pp['functions'][i]),
               pp['Ring']);
    return eq, B;

end proc:
EG := proc(ExplicitMatrix, n, m, rdim, x)
local k, v, c, l_c, l_c1, R, den;

    if convert(ExplicitMatrix[-1], 'set') = {0} then
        ExplicitMatrix[1..-2] :=
                            procname(ExplicitMatrix[1..-2], n, m, rdim-1, x);
        return ExplicitMatrix;
    end if;

    for k to rdim  do
        ExplicitMatrix[k] := Shift_Row(ExplicitMatrix[k], m, x);
    end do;

    v := LinearAlgebra:-NullSpace(ExplicitMatrix[1..-1, 1..m]^%T);
    if v = {} then
        return ExplicitMatrix;
    end if;

    l_c := 0;
    for k to rdim do
        l_c1 := Row_Length(ExplicitMatrix[k]);
        if v[1][k] <> 0 and l_c < l_c1 then
            c := k;
            l_c := l_c1;
        end if;
    end do;

    R := map(radnormal,(v[1]^%T).ExplicitMatrix);
    den := lcm(map(denom,convert(R, 'set'))[]);
    R := map(radnormal, den*R);
    if convert(R, 'set') = {0} then

        if c < rdim then
            ExplicitMatrix[c] := ExplicitMatrix[-1];
        end if;

        ExplicitMatrix[-1] := R;
        ExplicitMatrix[1..-2] :=
                            procname(ExplicitMatrix[1..-2], n, m, rdim-1, x);
        return ExplicitMatrix;
    else
        ExplicitMatrix[c] := R;
        return procname(ExplicitMatrix, n, m, rdim, x);
    end if;
end proc:

Shift_Row := proc(C, m, x)
local l, c, den;

    c := C;
    l := Row_Length(c);
    while l <> 0 and convert(c[1..m], 'set') = {0} do
        den := c[l];
        c := map(radnormal, c/den);
        c := <<c[m+1..-1] | Vector[row](m)>>[1] +
               map(diff, c, x);
        c := map(radnormal, c);
        den := lcm(map(denom,convert(c, 'set'))[]);
        c := map(radnormal, den*c);
        l := Row_Length(c);
    end do;

    return c;

end proc:
Row_Length := proc(c)
local k;
    for k from LinearAlgebra:-ColumnDimension(c) by -1 to 1 do
        if c[k] <> 0 then
            return k;
        end if;
    end do;
    return 0;
end proc:

check_arguments := proc()
local pp, m, m1, y, opts, S, M, vars, cfs, i, k, j, x, s_i, n, inverseA_l, t,
      dummy, dot_indts, indts, terms;
option remember;

    if nargs = 0 then
        error "The first argument (a system) is expected";
    end if;

    if args[1]::'Matrix' then
        if not args[1]::'Matrix'(square) then
            error "A square matrix is expected"
        end if;

        m := LinearAlgebra:-RowDimension(args[1]);
        if m = 0 then
            error "Non empty matrix is expected"
        end if;

        if nargs = 1 then
            error "The second argument (an independent variable) is expected"
        elif not args[2]::'name' then
            error "The second argument (a name) is expected"
        end if;
        x := args[2];

        if m = 1 then
            vars := [y(x)];
        else
            vars := [seq(y[i](x), i = 1..m)];
        end if;

        M := map(radnormal, args[1]);
        pp := table({
             'embracing_system' = <<Matrix(m, 'shape' = 'identity') | -M>>,
             'normal_form' = Array(0..0, [M]),
             'explicit_form' = Array(0..1, [-M, Matrix(m, 'shape' = 'identity')]),
             'order' = 1,
             'var' = x,
             'functions' = vars,
             'Ring' = OreTools:-SetOreRing(x, 'differential')});
        cfs := convert(args[1], 'set');

        opts := [args[3..-1]];

    elif args[1]::{'set', 'list'}({'algebraic', '`=`'}) then

        S := map(el->`if`(el::`=`, (lhs-rhs)(el), el),
                 convert(args[1], 'set')) minus {0};
        S := [S[]];
        if nops(S) = 0 then
            error "Non empty system is expected"
        end if;

        if nargs = 1 then
            error "The second argument (a list of unknowns) is expected"
        elif not args[2]::'list'('function'('name')) then
            error "The second argument (a list of function(name)) is expected"
        end if;
        vars := [args[2][]];
        y := map2(op, 0, vars);
        x := op(args[2][1]);
        m := nops(vars);
        m1 := nops(S);
        if nops({vars[]}) < m then
            error
           "The second argument (a list of different unknowns) is expected";
        end if;

        if {map(op, vars)[]} <> {x} then
            error
                 "The second argument (a list(function(%1))) is expected",
                                                                              x;
        end if;

        indts := select(el->has(x, el), y);
        if indts <> []  then
             error "%1 cannot be used for an unknown", indts[1](x);
        end if;

        pp := table({
             'var' = x,
             'functions' = vars,
             'Ring' = OreTools:-SetOreRing(x, 'differential')});

        opts := [args[3..-1]];

        if {eval(S, map(`=`, y, 0))[]} <> {0} then
            error "A homogeneous system is expected";
        end if;
        M := [seq(eval(S, map(`=`,subsop(i = NULL, y), 0)), i = 1 .. m)];
        try
            M := [seq(map(el ->
                     OreTools:-Converters:-FromLinearEquationToOrePoly(
                                           el+y[i](x)*dummy, y[i], pp['Ring']),
                          M[i]), i = 1 .. m)];
        catch:
            M := FAIL;
        end try;
        if has(M, FAIL) then
            error "A linear differential system is expected";
        end if;
        M := eval(M, dummy = 0);

        n := max(map(nops,map(op, M))) - 1;
        if n = 0 then
            error "A nonzero-order system is expected"
        end if;
        pp['order'] := n;

        if m > nops(S) then
            error "A full rank system is expected";
        end if;

        M := Matrix(nops(M[1]),m*(n+1),
                    {seq(seq(seq((k,(n+1-j+1)*m-(m-i))=op(j,M[i][k]),
                                 j = 1..nops(M[i][k])),
                             k = 1 .. nops(M[i])),
                         i = 1..m)});
        M := map(radnormal, M);
        if radnormal(expand({seq(S[k]-add(add(M[k][m*j-m+i]*
                  OreTools:-Apply(OrePoly(0$(n-j+1),1), vars[i], pp['Ring']),
                             j = 1..n+1),
                         i = 1..m),
                k = 1 .. nops(S))})) <> {0} then
            error "A linear system of %1 type is expected", op(2,pp['Ring']);
        end if;
        pp['explicit_form'] := Array(0..n,
                    [seq(M[1..-1,(n+1)*m-(m*(k+1)-1)..(n+1)*m-m*k], k = 0..n)]);
        cfs := convert(M, 'set');

    elif args[1]::{'algebraic',`=`} then
        try
            S := `if`(args[1]::`=`, (lhs-rhs)(args[1]), args[1]);
        catch:
            error "Wrong a system's representation";
        end try;

        if S::'Matrix(square)' and LinearAlgebra:-RowDimension(S) = 0 then
            error "Non empty system is expected"
        end if;

        if nargs = 1 then
            error "The second argument (an unknown) is expected"
        elif not args[2]::'function'('name') then
            error "The second argument (a function(name)) is expected"
        end if;
        y := op(0, args[2]);
        x := op(args[2]);
        if has(x, y) then
            error "%1 cannot be used for an unknown", y(x);
        end if;

        pp := table({
             'var' = x,
             'Ring' = OreTools:-SetOreRing(x, 'differential')});

        opts := [args[3..-1]];

        try
            M := convert(eval(S, y = 0), 'set');
        catch:
            error "Wrong a system's representation";
        end try;
        if M = {} then
            error "Non empty matrices are expected"
        elif M <> {0} then
            error "A homogeneous linear system with matrix coefficients is expected";
        end if;

        dot_indts := indets(S, `.`);
        dot_indts, indts := selectremove(el -> nops(el) = 2 and
                                               op(1,el)::'Matrix', dot_indts);
        if indts <> {} then
            error "Wrong a system's representation";
        end if;

        M := eval(S, {seq(dot_indts[i] = cfs[i]*op(2,dot_indts[i]),
                                                    i = 1 .. nops(dot_indts))});
        try
            M := OreTools:-Converters:-FromLinearEquationToOrePoly(
                                              M + dummy*y(x), y, pp['Ring']);
        catch:
            M := FAIL;
        end try;
        if has(M, FAIL) then
            error "A linear system of %1 type is expected", op(2,pp['Ring']);
        end if;
        M := eval(M, dummy = 0);
        n := nops(M) - 1;
        if n = 0 then
            error "A nonzero-order system is expected"
        end if;

        m := map(el-> [LinearAlgebra:-Dimension(op(1, el))], dot_indts);
        if nops(m) <> 1 then
            error "Wrong a system's representation";
        elif m[1][1] < m[1][2] then
            error "A full rank system is expected";
        end if;
        m1, m := m[1][];
        if m = 0 then
            error "Non empty matrices are expected"
        elif m = 1 then
            vars := [y(x)];
        else
            vars := [seq(y[i](x), i = 1..m)];
        end if;
        pp['functions'] := vars;
        pp['order'] := n;

        try
            dot_indts := map(el ->
               [op(1, el),
                OreTools:-Converters:-FromLinearEquationToOrePoly(
                                                     op(2, el), y, pp['Ring']),
                coeff(S,el)],
                             dot_indts);
        catch:
            error "Wrong a system's representation";
        end try;

        if select(el -> nops(el[2]) = 0 or
                        op(-1, el[2]) <> 1 or
                        not {op(1..-2, el[2])} in {{},{0}}, dot_indts) <> {} then
            error "Wrong a system's representation";
        end if;
        M := Matrix(m1,m*(n+1));
        for i in dot_indts do
            M[1..-1,m*(n+1)-m*nops(i[2])+1..m*(n+1)-m*nops(i[2])+m] :=
                 M[1..-1,m*(n+1)-m*nops(i[2])+1..m*(n+1)-m*nops(i[2])+m] +
                 i[1]*i[3];
        end do;
        M := map(radnormal, M);
        pp['explicit_form'] := Array(0..n,
                    [seq(M[1..-1,(n+1)*m-(m*(k+1)-1)..(n+1)*m-m*k], k = 0..n)]);
        cfs := convert(M, 'set');

    else
        error
         "The first argument (a system in matrix or list/set form) is expected";
    end if;

    if has(cfs, y) then
        if y::'name' then
            error "%1 cannot be used for an unknown", y(x);
        else
           y := select(el->has(cfs, el), y);
           if nops(y) = 1 then
                error "The name %1 cannot be used for unknowns", y[];
           else
                error "The names %1 cannot be used for unknowns", y;
           end if;
        end if;
    end if;

    # check system's coefficients
    if not cfs::'set'('ratpoly'('anything',x)) then
        error
         "The system's coefficients of type ratpoly(anything, %1) are expected",
                                                                              x;
    end if;

    if not assigned(pp['embracing_system']) then
        M := EG(M, n, m, m1, pp['var']);

        if m < m1 and convert(M[m+1..-1], 'set') <> {0} then
            error "It's impossible"
        end if;
        if convert(M[m], 'set') = {0} then
            error "A full rank system is expected";
        end if;
        pp['embracing_system'] := M[1..m];
        inverseA_l := (M[1 .. m,1 .. m])^(-1);
        pp['normal_form'] := Array(0..n-1,
            [seq(copy(map(radnormal,
                   -inverseA_l.M[1 .. m, (n+1)*m-(m*(k+1)-1)..(n+1)*m-m*k])),
                 k=0..n-1)]);
    end if;

    if nops(opts) > 0 and opts[1]::'name' then
        t := opts[1];
        opts := opts[2..-1];
    else
        t := FAIL;
    end if;
    if not opts::'list'(`=`) then
        error
    "optional arguments are expected to be of type equation, but received '%1'",
                                                      remove(type, opts, `=`)[];
    end if;

    terms := 0;
    hasoption(opts, 'order'='nonnegint', 'terms', 'opts');

    s_i := `if`(args[2]::'list', vars[1], 1);
    hasoption(opts, 'select_indicator', 's_i', 'opts');
    if not s_i::'posint' and not s_i::'function'('name') then
        i := s_i(args[1], args[2]);
    else
        i := s_i;
    end if;
    if i::'posint' then
        if args[2]::'set' then
            error "'select_indicator' is expected a function from %1", vars;
        elif i > m then
            error "The option argument 'select_indicator' is expected <= %1", m;
        end if;
    elif i::'function'('name') then
        if not member(i, vars, 'i') then
            error "'select_indicator' is expected a function from %1", vars;
        end if;
    end if;

    if args[1]::'Matrix' then
        input_type := 'Matrix';
    else
        input_type := 'list';
    end if;

    if opts <> [] then
        error "Invalid optional arguments; first one is '%1'", opts[1];
    end if;

    return pp, t, terms, i, s_i;

end proc:

    try ModuleApply(args)
    catch:
        error;
    end try;

end proc: