ode2singular.mac

/*Andreas Hermann 2020*/
/*The function ode2singular tries to find the general solution y(x) to the equation p2(x) * y''(x) + p1(x) * y'(x) + p0(x) * y(x) = 0.

Example of usage:
p2: x^2$
p1: x$
p0: 1$
eq: p2 * 'diff(y, x, 2) + p1 * 'diff(y, x) + p0 * y = 0$
ode2singular(eq, y, x);

In order to solve the equation with a change of variable replacing x by phi(x):
p2: 1$
p1: 0$
p0: a - 2 * q * cos(2 * x)$
phi: cos(x)$
phi_inv: acos(x)$
eq: p2 * 'diff(y,x,2) + p1 * 'diff(y,x) + p0 * y = 0$
ode2singular_transform(eq, y, x, phi, phi_inv);

We write y(x) = h(x) * g(x) such that g(x) satsfies an equation of the form g''(x) + q0(x) * g(x) = 0.
We define args_list: [p1 / p2, p0 / p2, q0, argument, inhom]
where argument is the current variable and inhom is the inhomogeneity of the equation.
At the beginning, we have argument = x.
At the moment, some inhomogeneous equations can be solved.
We use some equations from the Digital Library of Mathematical Functions
https://dlmf.nist.gov/
and we put references as we use them.
*/


ode2singular(eq, dep_var, x) := block(
    [radexpand: all, eq_new, p2, p1, p0, inhom, h_function, q0, argument, args_list, answer],
    eq_new: expand(lhs(eq) - rhs(eq)),
    p2: coeff(eq_new, 'diff(dep_var, x, 2)),
    p1: coeff(eq_new, 'diff(dep_var, x)),
    p0: coeff(eq_new, dep_var),
    inhom: ratsimp(eq_new - p2 * 'diff(dep_var, x, 2) - p1 * 'diff(dep_var, x) - p0 * dep_var),
    if freeof(dep_var, inhom) = false or ratsimp(p2) = 0 then
        (
        print("The function ode2singular cannot solve this equation since it is not linear of order 2"),
        /* homogeneous */
        return([])
        )
    else
        (
        p1: ratsimp(p1 / p2),
        p0: ratsimp(p0 / p2),
        inhom: ratsimp(-inhom / p2),
        h_function: radcan(exp(-integrate(p1, x) / 2)),
        q0: trigsimp(ratsimp(-diff(p1, x) / 2 - p1^2 / 4 + p0)),
        /*print("This function tries to find the general solution f to the equation f'' + p1 * f' + p0 * f = 0"),
        print("p1 = ", factor(p1)),
        print("p0 = ", factor(p0)),
        print("The dependent variable is", dep_var, ", the independent variable is", x),
        print("Every solution has the form f = h * g, where h = ", h_function),
        print("and g is a solution to g'' + Q * g = 0, where"),
        print("Q = ", factor(q0)),*/
        argument: x,
        args_list: [p1, p0, q0, argument, inhom],
        answer: determine_singular_points(args_list, h_function, dep_var, x),
        return(answer)
        )
    )$


ode2singular_transform(eq, dep_var, x, phi, phi_inv) := block(
    [radexpand: all, eq_new, p2, p1, p0, inhom, h_function, q0, argument, args_list, args_new, answer],
    eq_new: expand(lhs(eq) - rhs(eq)),
    p2: coeff(eq_new, 'diff(dep_var, x, 2)),
    p1: coeff(eq_new, 'diff(dep_var, x)),
    p0: coeff(eq_new, dep_var),
    inhom: ratsimp(eq_new - p2 * 'diff(dep_var, x, 2) - p1 * 'diff(dep_var, x) - p0 * dep_var),
    if freeof(dep_var, inhom) = false or ratsimp(p2) = 0 then
        (
        print("The function ode2singular cannot solve this equation since it is not linear of order 2"),
        return([])
        )
    else
        (
        p1: ratsimp(p1 / p2),
        p0: ratsimp(p0 / p2),
        inhom: ratsimp(-inhom / p2),
        h_function: radcan(exp(-integrate(p1, x) / 2)),
        q0: trigsimp(ratsimp(-diff(p1, x) / 2 - p1^2 / 4 + p0)),
        /*print("This function tries to find the general solution f to the equation f'' + p1 * f' + p0 * f = 0"),
        print("p1= ", factor(p1)),
        print("p0= ", factor(p0)),
        print("The dependent variable is", dep_var, ", the independent variable is", x),
        print("Every solution has the form f = h * g, where h = ", h_function),
        print("and g is a solution to g'' + Q * g = 0, where"),
        print("Q = ", factor(q0)),*/
        argument: x,
        args_list: [p1, p0, q0, argument, inhom],
        args_new: inv_transform(args_list, phi, phi_inv, x),
        args_new[2]: ratsimp(trigexpand(args_new[2])),
        args_new[1]: ratsimp(trigexpand(args_new[1])),
        args_new[3]: ratsimp(trigexpand(args_new[3])),
        args_new[4]: ratsimp(trigexpand(args_new[4])),
        /* maybe include the lines "args_new[i]: ratsimp(trigexpand(args_new[i]))," in the function inv_transform ??? */
        /*print("new list: ", args_new),*/
        answer: determine_singular_points(args_new, h_function, dep_var, x),
        return(answer)
        )
    )$



determine_singular_points(args_list,h_function,dep_var,x) := block(
    [i, fact_data, reg_sing, irreg_sing, reg_mult, irreg_mult, lim_inf, lim_minf, answer],
    /*Replaced denom by ratdenom in the following line since some denominators were not determined correctly*/
    fact_data: determine_factorization(ratdenom(args_list[3]), x),
    /*We have fact_data = [coeff, sing_points, multiplicities]*/
    if fact_data[1]#inf then
        (
        reg_sing: [],
        irreg_sing: [],
        reg_mult: [],
        irreg_mult: [],
        /*-----------------------------------------------------------------------------------------------
        decide for every singular point of the equation in Q-form whether it is regular or 
        irregular: If q0 has a singularity of order at most 2 then the singular point is regular, 
        otherwise it is irregular.
        -----------------------------------------------------------------------------------------------*/
        for i:1 thru length(fact_data[2]) do
            (
            if fact_data[3][i] > 2 then
                (
                push(fact_data[2][i], irreg_sing),
                push(fact_data[3][i], irreg_mult)
                )
            else
                (
                push(fact_data[2][i], reg_sing),
                push(fact_data[3][i], reg_mult)
                )
            ),
        /*-----------------------------------------------------------------------------------------------
        Infinity is always a singular point of the equation in Q-form.
        Decide whether infinity is a regular or an irregular singular point.
        If x^2*q0(x) is bounded as abs(x) goes to infinity then infinity is a regular singular 
        point, otherwise it is an irregular singular point.
        x=inf will be the first item in the list reg_sing or irreg_sing respectively
        -----------------------------------------------------------------------------------------------*/
        lim_inf: abs(limit(x^2*args_list[3], x, inf)),
        lim_minf: abs(limit(x^2*args_list[3], x, minf)),
        if testmode=true then
            (
            print("q0=", args_list[3]),
            print("lim_inf=", lim_inf),
            print("lim_minf=", lim_minf)
            ),
        if lim_inf#inf and lim_minf#inf and lim_inf#und and lim_minf#und then
            (
            push(x = inf, reg_sing),
            /*We don't determine the multiplicity of infinity.*/
            push(unknown, reg_mult)
            )
        else
            (
            push(x = inf, irreg_sing),
            /*We don't determine the multiplicity of infinity.*/
            push(unknown, irreg_mult)
            ),
        if length(reg_sing) = 0 then
            (
            print("The equation has no regular singular points.")
            )
        else
            (
            print("The equation has the regular singular points ", reg_sing)
            /*print("with multiplicities ", reg_mult)*/
            ),
        if length(irreg_sing) = 0 then
            (
            print("The equation has no irregular singular points.")
            )
        else
            (
            print("The equation has the irregular singular points ", irreg_sing)
            /*print("with multiplicities ", irreg_mult)*/
            ),
        /* constant coefficients */
        if diff(args_list[3], x) = 0 then 
            (
            answer: solve_const_coeff(args_list, h_function, dep_var, x)
            )
        /* 1 singular point */
        else if length(reg_sing) = 0 and length(irreg_sing) = 1 then
            (
            answer: solve_parabolic_cylinder(args_list, h_function, dep_var, x)
            )
        /* 2 singular points */
        else if length(reg_sing) = 1 and length(irreg_sing) = 1 then
            (
            answer: solve_bessel_or_kummer(args_list, h_function, dep_var, x, reg_sing, irreg_sing)
            )
        else if length(reg_sing) = 2 and length(irreg_sing) = 0 then
            (
            answer: solve_power(args_list, h_function, dep_var, x, reg_sing, irreg_sing)
            )
        /* 3 singular points */
        else if length(reg_sing) = 3 and length(irreg_sing) = 0 then 
            (
            answer: solve_legendre_or_hypergeo(args_list, h_function, dep_var, x, reg_sing, irreg_sing, reg_mult)
            )
        else if length(reg_sing) = 2 and length(irreg_sing) = 1 then
            (
            answer: solve_bessel_or_kummer_mobius_or_spheroid(args_list, h_function, dep_var, x, reg_sing, irreg_sing)
            )
        /* 4 singular points */
        else if length(reg_sing) = 4 and length(irreg_sing) = 0 then
            (
            answer: solve_4_reg_sing(args_list, h_function, dep_var, x, reg_sing, irreg_sing)
            )
        else
            (
            print("Equations with that number of singular points have not been implemented so far."),
            answer: []
            )
        )
    else
        (
        print("Error: Could not factor the denominator of Q=",factor(args_list[3])),
        answer: []
        ),
    if answer=[] then
        (
        print("The function ode2singular could not solve the equation."),
        print("You can try to transform the equation to a new independent variable new_var by calling"),
        print("ode2singular_transform(eq, y, x, new_var, new_var_inv) where"),
        print("eq is your equation, y is the dependent variable (in your case", dep_var, "),"),
        print("x is the independent variable (in your case", x,")"),
        print("new_var is the new variable written as a function of ", x),
        print("and new_var_inv is the inverse function of new_var, also written as a function of ", x)
        /*print("new_var and new_var_inv should both be written as functions of", x),
        print("and new_var_inv is the inverse function of new_var.")*/
        ),
    return(answer)
    )$


complex_sqrt(re, im) := block([sol_re, sol_im],
    /*This function computes a complex square root of the number re + %i * im*/
    sol_re: radcan(sqrt(ratsimp((radcan(sqrt(re^2+im^2))+re)/2))),
    sol_im: radcan(sqrt(ratsimp((radcan(sqrt(re^2+im^2))-re)/2))),
    if im<0 then
        (
        sol_im: -sol_im
        ),
    return(sol_re+%i*sol_im)
    )$


simpl(f) := block([],
    return(factor(ratsimp(radcan(f))))
    )$


/* remove multiplicative constants */
remove_constants(f,x) := block([logexpand: all, g],
    if f=0 then
        (
        return(0)
        )
    else
        (
        /*g: factor(ratsimp(radcan(f))),
        g: diff(log(g),x),*/
        g: diff(log(f),x),
        return(radcan(exp(integrate(g,x))))
        )
    )$


write_solution(args_list, dep_var, x, f_1, f_2, wronski, answer) := block([A, B, sol],
    if wronski = tbc and args_list[5] # 0 then
        (
        print("The equation is inhomogeneous. The solution to this equation has not yet been implemented."),
        /*print("The general solution to the homogeneous equation is:"),
        push(dep_var=ratsimp(trigsimp( %k1 * remove_constants(f_1, x) + %k2 * remove_constants(f_2, x))), answer),*/
        return(answer)
        )
    else
        (
        /* formula for variation of parameters, see e.g. https://en.wikipedia.org/wiki/Variation_of_parameters */
        A: -integrate(f_2*args_list[5]/wronski,x),
        B: integrate(f_1*args_list[5]/wronski,x),
        sol: A * f_1 + B * f_2 + %k1 * remove_constants(f_1, x) + %k2 * remove_constants(f_2, x),
        /* try simplifications  */
        /* sol: ratsimp(sol), */
        push(dep_var=sol, answer),
        return(answer)
        )
    )$


solve_const_coeff(args_list, h_function, dep_var, x) := block(
    [re, im, sol, h_function_const, m_function, f_1, f_2, wronski, answer, factor],
    answer: [],
    if diff(args_list[3],x)#0 then
        (
        print("Error: The function solve_const_coeff should only be called if Q is constant.")
        )
    else
        (
        re: realpart(args_list[3]),
        im: imagpart(args_list[3]),
        h_function_const: 1,
        m_function: diff(args_list[4], x)^(1/2) * subst(args_list[4], x, h_function_const),
        factor: h_function / m_function,
        if asksign(im)=zero then
            (
            if asksign(re)=pos then
                (
                f_1: factor * cos(sqrt(re)*args_list[4]),
                f_2: factor * sin(sqrt(re)*args_list[4]),
                wronski: factor^2 * sqrt(re)*diff(args_list[4],x)
                )
            else if asksign(re)=neg then 
                (
                f_1: factor * exp(sqrt(-re)*args_list[4]),
                f_2: factor * exp(-sqrt(-re)*args_list[4]),
                wronski: -2 * factor^2 * sqrt(-re) * diff(args_list[4],x)
                )
            else
                (
                f_1: factor * args_list[4],
                f_2: factor,
                wronski: -factor^2 * diff(args_list[4],x)
                )
            )
        else
            (
            sol: complex_sqrt(-re,-im),
            f_1: factor * exp(sol*args_list[4]),
            f_2: factor * exp(-sol*args_list[4]),
            wronski: -2 * factor^2 * sol*diff(args_list[4],x)
            ),
        answer: write_solution(args_list, dep_var, x, f_1, f_2, wronski, answer)
        ),
    return(answer)
    )$


solve_power(args_list, h_function, dep_var, x, reg_sing, irreg_sing) := block(
    [beta, coeffC, h_function_power, phi, phi_inv, args_new, m_function_new, discrim, f_1, f_2, wronski, answer, factor],
    answer: [],
    if length(reg_sing)#2 or length(irreg_sing)#0 then 
        (
        print("Error: The function solve_power should only be called in case of 2 regular singular points and 0 irregular singular points.")
        )
    else if rhs(reg_sing[1])#inf and rhs(reg_sing[2])#inf then
        (
        print("Error: The function solve_power should only be called if inf is a regular singular point.")
        )
    else
        (
        if rhs(reg_sing[1])=inf then 
            (
            beta: rhs(reg_sing[2])
            )
        else
            (
            beta: rhs(reg_sing[1])
            ),
        coeffC: ratsimp(args_list[3]*(x-beta)^2),
        if diff(coeffC,x)#0 then
            (
            print("Error: The function solve_power should not be called if Q(x) is not of the form C/(x-beta)^2 with C constant.")
            )
        else
            (
            /* The solutions are of the form x^k where k can be a complex number */
            h_function_power: 1,
            phi: x-beta,
            phi_inv: x+beta,
            args_new: inv_transform(args_list, phi, phi_inv, x),
            m_function_new: diff(args_new[4], x)^(1/2) * subst(args_new[4], x, h_function_power),
            factor: h_function / m_function_new,
            if coeffC=1/4 then
                (
                f_1: factor * sqrt(args_new[4]),
                f_2: factor * sqrt(args_new[4])*log(args_new[4]),
                wronski: factor^2 * diff(args_new[4], x)
                )
            else
                (
                discrim: complex_sqrt(1-4*realpart(coeffC),-4*imagpart(coeffC)),
                f_1: factor * args_new[4]^((1+realpart(discrim))/2)*(cos(imagpart(discrim)/2*log(args_new[4]))+%i*sin(imagpart(discrim)/2*log(args_new[4]))),
                f_2: factor * args_new[4]^((1-realpart(discrim))/2)*(cos(imagpart(discrim)/2*log(args_new[4]))-%i*sin(imagpart(discrim)/2*log(args_new[4]))),
                wronski: -discrim * factor^2 * diff(args_new[4], x)
                ),
            answer: write_solution(args_new, dep_var, x, f_1, f_2, wronski, answer)
            )
        ),
    return(answer)
    )$


solve_parabolic_cylinder(args_list, h_function, dep_var, x) := block(
    [answer, f_1, f_2, w0, wronski, det_data, coeffA, coeffC, exponent, x0, alpha, k, sol_a, sol_b, sol_b_2, lambda, phi, phi_inv, args_new, m_function_new, h_function_airy, h_function_bessel, h_function_hermite, h_function_kummer, factor],
    answer: [],
    /* check whether q0(x) is of the form A*(x-x0)^exponent   */
    det_data: detect_power_1(args_list[3],x),
    coeffA: det_data[1],
    exponent: det_data[2],
    x0: det_data[3],
    if exponent=1 then
        (
        /*-----------------------------------------------
        Airy functions Ai(alpha*(x-x0)), Bi(alpha*(x-x0)), q0(x)=-alpha^3*(x-x0)
        Need to solve alpha^3=-coeffA for alpha
        -----------------------------------------------*/
        alpha: (-coeffA)^(1/3),
        h_function_airy: 1,
        phi: alpha*(x-x0),
        phi_inv: x/alpha+x0,
        args_new: inv_transform(args_list, phi, phi_inv, x),
        m_function_new: diff(args_new[4],x)^(1/2) * subst(args_new[4], x, h_function_airy),
        print("The solutions can be expressed in terms of Airy functions."),
        factor: h_function / m_function_new,
        f_1: factor * airy_ai(args_new[4]),
        f_2: factor * airy_bi(args_new[4]),
        /*Wronskian of Airy functions see https://dlmf.nist.gov/9.2.E7  */
        w0: 1 / %pi,
        wronski: factor^2 * diff(args_new[4], x) * subst(args_new[4], x, w0),
        answer: write_solution(args_new, dep_var, x, f_1, f_2, wronski, answer)
        ),
    if exponent#inf and exponent#-2 then
        (
        /*--------------------------------------------------
        Bessel functions: q0(x)=alpha^2*k^2*(x-x0)^exponent.
        Equations to solve for alpha, sol_a, k:
        exponent=2*k-2,
        alpha^2*k^2=coeffA,
        1-4*sol_a^2*k^2=0
        --------------------------------------------------*/
        k: (exponent+2)/2,
        sol_a: 1/(2*k),
        alpha: sqrt(coeffA)/k,
        h_function_bessel: x^(-1/2),
        phi: alpha*(x-x0)^k,
        phi_inv: (x/alpha)^(1/k)+x0,
        args_new: inv_transform(args_list, phi, phi_inv, x),
        m_function_new: diff(args_new[4], x)^(1/2) * subst(args_new[4], x, h_function_bessel),
        answer: bessel_fund_system(args_new, h_function/m_function_new, dep_var, x, sol_a, answer)
        ),
    if exponent#-2 then
        (
        /* check whether q0(x) is of the form  coeffA*(x-x0)^(2*exponent-2)+coeffC*(x-x0)^(exponent-2) */
        det_data: detect_power_4(args_list[3], x),
        coeffA: det_data[1],
        coeffC: det_data[2],
        exponent: det_data[3],
        x0: det_data[4],
        if exponent#inf and exponent#0 and coeffA#0 then
            (
            /*-------------------------------------------------------------------
            Hermite equation: H_lambda(alpha*(x-x0)), lambda non-negative integer
            q0(x)=-alpha^4*(x-x0)^2+alpha^2*(1+2*lambda)
            -------------------------------------------------------------------*/
            lambda: ratsimp(radcan((coeffC/sqrt(-coeffA)-1)/2)),
            if exponent=2 and testmode#true and askinteger(lambda)=yes and asksign(lambda)#neg then
                (
                print("The solutions can be expressed in terms of Hermite polynomials"),
                print("of degree ", lambda, " and confluent hypergeometric functions."),
                alpha: (-coeffA)^(1/4),
                h_function_hermite: exp(x^2/2),
                phi: alpha*(x-x0),
                phi_inv: x/alpha+x0,
                args_new: inv_transform(args_list, phi, phi_inv, x),
                m_function_new: diff(args_new[4],x)^(1/2)*subst(args_new[4],x,h_function_hermite),
                factor: h_function / m_function_new,
                if askinteger(lambda,even)=yes then
                    (
                    f_1: factor * args_new[4] * kummer_m((1-lambda)/2, 3/2, args_new[4]^2),
                    f_2: factor * hermite(lambda,args_new[4]),
                    /* hermite(lambda,x) is equal to 2^lambda*x*kummer_u((1-lambda)/2, 3/2, x^2) by https://dlmf.nist.gov/18.11.E3   */
                    /* https://dlmf.nist.gov/13.2.E34 */
                    w0: -x^(-3/2) * exp(x) * gamma(3/2) / gamma((1-lambda)/2),
                    wronski: (factor * args_new[4])^2 * 2^lambda * diff(args_new[4]^2, x) * subst(args_new[4]^2, x, w0)
                    )
                else
                    (
                    f_1: factor * kummer_m(-lambda/2,1/2,args_new[4]^2),
                    f_2: factor * hermite(lambda,args_new[4]),
                    /* hermite(lambda,x) is equal to 2^lambda*kummer_u(-lambda/2,1/2,x^2) by https://dlmf.nist.gov/18.11.E3   */
                    /* https://dlmf.nist.gov/13.2.E34 */
                    w0: -x^(-1/2) * exp(x) * gamma(1/2) / gamma(-lambda/2),
                    wronski: factor^2 * 2^lambda * diff(args_new[4]^2, x) * subst(args_new[4]^2, x, w0)
                    ),
                answer: write_solution(args_new, dep_var, x, f_1, f_2, wronski, answer)
                )
            else
                (
                /*------------------------------------------------------
                Kummer equation: q0(x)=-alpha^2*exponent^2/4*(x-x0)^(2*exponent-2)+(2*sol_b-4*sol_a)/4*alpha*exponent^2*(x-x0)^(exponent-2).
                Equations to solve for alpha, sol_a, sol_b:
                -alpha^2*exponent^2/4=coeffA
                (2*sol_b-4*sol_a)*alpha*exponent^2/4=coeffC
                1-(sol_b-1)^2*exponent^2=0
                The last equation has 2 solutions sol_b, sol_b_2.
                -------------------------------------------------------*/
                sol_b: ratsimp(1+1/exponent),
                sol_b_2: ratsimp(1-1/exponent),
                /* first consider sol_b */
                alpha: ratsimp(2*radcan(sqrt(-coeffA))/exponent),
                sol_a: ratsimp(sol_b/2-coeffC/(alpha*exponent^2)),
                h_function_kummer: x^(-sol_b/2)*exp(x/2),
                phi: alpha*(x-x0)^exponent,
                phi_inv: (x/alpha)^(1/exponent)+x0,
                args_new: inv_transform(args_list, phi, phi_inv, x),
                m_function_new: diff(args_new[4],x)^(1/2)*subst(args_new[4],x,h_function_kummer),
                print("The solutions can be expressed in terms of parabolic cylinder functions."),
                answer: kummer_fund_system(args_new, h_function/m_function_new, dep_var, x, sol_a, sol_b, answer),
                /* consider also the solution sol_b_2 */
                sol_a: ratsimp(sol_b_2/2-coeffC/(alpha*exponent^2)),
                h_function_kummer: x^(-sol_b_2/2)*exp(x/2),
                m_function_new: diff(args_new[4],x)^(1/2)*subst(args_new[4],x,h_function_kummer),
                answer: kummer_fund_system(args_new, h_function/m_function_new, dep_var, x, sol_a, sol_b_2, answer)
                )
            )
        else
            (
            print("Error: The function_parabolic_cylinder cannot find a solution.")
            )
        ),
    return(answer)
    )$


bessel_fund_system(args_list, factor, dep_var, x, sol_a, answer) := block(
    [f_1, f_2, w0, wronski],
    print("The solutions can be expressed in terms of Bessel functions"),
    if testmode=true then
        (
        f_1: factor * bessel_j(sol_a, args_list[4]),
        f_2: factor * bessel_y(sol_a, args_list[4]),
        wronski: tbc
        )
    else
        (
        f_1: factor * bessel_j(sol_a, args_list[4]),
        f_2: factor * bessel_y(sol_a, args_list[4]),
        /* https://dlmf.nist.gov/10.5.E2  */
        w0: 2/(%pi*x),
        wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
        ),
    answer: write_solution(args_list, dep_var, x, f_1, f_2, wronski, answer),
    return(answer)
    )$


/*The function kummer_fund_system assigns
dep_var=%k1*factor*kummer_1(sol_a,sol_b,args_list[4])+%k2*factor*kummer_2(sol_a,sol_b,args_list[4])
where kummer_1 and kummer_2 are two linearly independent solutions of the Kummer equation.
The functions kummer_1 and kummer_2 are the functions kummer_m or kummer_u or related functions or generalized Laguerre polynomials.*/

kummer_fund_system(args_list, factor, dep_var, x, sol_a, sol_b, answer) := block(
    [expand_hypergeometric: true, f_1, f_2, wronski, beta],
    if askinteger(sol_a)=yes and asksign(sol_a)#pos then
        (
        print("The solutions can be expressed in terms of generalized Laguerre polynomials"),
        print("of degree ", -sol_a, " with parameter ", ratsimp(sol_b-1), " and confluent hypergeometric functions."),
        if askinteger(sol_b)=no then 
            (
            f_1: factor * gen_laguerre(-sol_a, ratsimp(sol_b-1), args_list[4]),
            f_2: factor * args_list[4]^(ratsimp(1-sol_b)) * kummer_m(ratsimp(1+sol_a-sol_b), ratsimp(2-sol_b), args_list[4]),
            /* https://dlmf.nist.gov/13.2.E33 */
            w0: sin(%pi*sol_b) * x^(-sol_b) * exp(x) / %pi * gamma(sol_b) * gamma(2-sol_b),
            /* gen_laguerre(-a,b-1,x) is equal to pochhammer(b,-a)/(-a)! * kummer_m(a,b,x) by   https://dlmf.nist.gov/18.11.E2 */
            wronski: pochhammer(sol_b,-sol_a)/(-sol_a)! * factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )
        /* the function gen_laguerre(n, alpha, x) returns a division by 0 error e.g. for n=1, alpha=-1  */
        else
            (
            f_1: factor * subst(ratsimp(sol_b-1), beta, ratsimp(gen_laguerre(-sol_a, beta, args_list[4]))),
            f_2: factor * exp(args_list[4]) * kummer_u(ratsimp(sol_b-sol_a), sol_b, -args_list[4]),
            /*   https://dlmf.nist.gov/13.2.E38  */
            w0: (-1)^(sol_a-sol_b) * x^(-sol_b) * exp(x),
            /* gen_laguerre(-a,b-1,x) is equal to (-1)^(-a)/(-a)! * kummer_u(a,b,x) by   https://dlmf.nist.gov/18.11.E2 */
            wronski: (-1)^(-sol_a)/(-sol_a)! * factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            ),
        answer: write_solution(args_list, dep_var, x, f_1, f_2, wronski, answer)
        ),
    if testmode=true then
        (
        print("The solutions can be expressed in terms of confluent hypergeometric functions."),
        f_1: factor * kummer_m(sol_a, sol_b, args_list[4]),
        f_2: factor * kummer_u(sol_a, sol_b, args_list[4]),
        wronski: tbc
        )
    else 
        (
        print("The solutions can be expressed in terms of confluent hypergeometric functions."),
        /*  sol_b in {1, 2, 3, ... }  */
        if askinteger(sol_b)=yes and asksign(sol_b)=pos then
            (
            f_1: factor * hypergeometric([sol_a], [sol_b], args_list[4]),
            /*  sol_a in {1, 2, ..., sol_b-1}  */
            if askinteger(sol_a)=yes and asksign(sol_a)=pos and asksign(sol_a-sol_b)=neg then
                (
                /*  https://dlmf.nist.gov/13.2.E29  */
                f_2: sum((k-1)!*args_list[4]^(-k)/((sol_b-1-k)!*(k-sol_a)!), k, sol_a, sol_b-1),
                f_2: factor * f_2,
                /*  we have diff(f_2, x) = -(sol_b-1)!/(sol_b-1-sol_a)! * x^(-sol_b) + o(x^(-sol_b)) as x->0  */
                w0: -(sol_b-1)! / (sol_b-1-sol_a)! * x^(-sol_b) * exp(x),
                wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
                )
            /*  sol_a in {0, -1, -2, ...} */
            else if askinteger(sol_a)=yes and asksign(sol_a)#pos then
                (
                /*  https://dlmf.nist.gov/13.2.E28  */
                f_2: sum((sol_b-1)!*(k-1)!*args_list[4]^(-k)/((sol_b-1-k)!*pochhammer(1-sol_a,k)), k, 1, sol_b-1)
                    -sum(pochhammer(sol_a,k)*args_list[4]^k/(pochhammer(sol_b,k)*k!)*(log(args_list[4])+psi[0](1-sol_a-k)-psi[0](1+k)-psi[0](sol_b+k)), k, 0, -sol_a)
                    +(-1)^(1-sol_a)*(-sol_a)!*sum((k-1+sol_a)!*args_list[4]^k/(pochhammer(sol_b,k)*k!), k, 1-sol_a, inf),
                f_2: factor * f_2,
                /*  we have diff(f_2, x) = -(sol_b-1)! * (sol_b-1)! / pochhammer(1-sol_a,sol_b-1) * x^(-sol_b) + o(x^(-sol_b)) as x->0  */
                w0: -(sol_b-1)! * (sol_b-1)! / pochhammer(1-sol_a,sol_b-1) * x^(-sol_b) * exp(x),
                wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
                )
            else
                (
                /*  https://dlmf.nist.gov/13.2.E27  */
                f_2: sum((sol_b-1)!*(k-1)!*args_list[4]^(-k)/((sol_b-1-k)!*pochhammer(1-sol_a,k)), k, 1, sol_b-1)
                    -sum(pochhammer(sol_a,k)*args_list[4]^k/(pochhammer(sol_b,k)*k!)*(log(args_list[4])+psi[0](sol_a+k)-psi[0](1+k)-psi[0](sol_b+k)), k, 0, inf),
                f_2: factor * f_2,
                /*  we have diff(f_2, x) = -(sol_b-1)! * (sol_b-1)! / pochhammer(1-sol_a,sol_b-1) * x^(-sol_b) + o(x^(-sol_b)) as x->0  */
                w0: -(sol_b-1)! * (sol_b-1)! / pochhammer(1-sol_a,sol_b-1) * x^(-sol_b) * exp(x),
                wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
                )
            )
        /*  sol_b in {0, -1, -2, ...}  */
        else if askinteger(sol_b)=yes and asksign(sol_b)#pos then
            (
            f_1: factor * args_list[4]^(1-sol_b) * hypergeometric([ratsimp(1+sol_a-sol_b)], [ratsimp(2-sol_b)], args_list[4]),
            /*  sol_a in {0, -1, -2, ..., sol_b}  */
            if askinteger(sol_a)=yes and asksign(sol_a)#pos and asksign(sol_a-sol_b)#neg then
                (
                /*  https://dlmf.nist.gov/13.2.E32  */
                f_2: sum((k-1)!*args_list[4]^(1-sol_b-k)/((1-sol_b-k)!*(k+sol_b-sol_a-1)!), k, sol_a-sol_b+1, 1-sol_b),
                f_2: factor * f_2,
                /*  we have -f_2 * diff(f_1, x) = -(1-sol_b)! / (-sol_a)! * x^(-sol_b) + o(x^(-sol_b)) as x->0  */
                w0: -(1-sol_b)! / (-sol_a)! * x^(-sol_b) * exp(x),
                wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
                )
            /*  sol_a in {sol_b-1, sol_b-2, ...}  */
            else if askinteger(sol_a)=yes and asksign(sol_a-sol_b)=neg then
                (
                /*  https://dlmf.nist.gov/13.2.E31  */
                f_2: sum((1-sol_b)!*(k-1)!*args_list[4]^(1-sol_b-k)/((1-sol_b-k)!*pochhammer(sol_b-sol_a,k)), k, 1, 1-sol_b)
                    -sum(pochhammer(sol_a-sol_b+1,k)*args_list[4]^(1-sol_b+k)/(pochhammer(2-sol_b,k)*k!)
                        *(log(args_list[4])+psi[0](sol_b-sol_a-k)-psi[0](1+k)-psi[0](2-sol_b+k)), k, 0, sol_b-sol_a-1)
                    +(-1)^(-sol_b-sol_a)*(sol_b-sol_a-1)!*sum((k+sol_a-sol_b)!/(pochhammer(2-sol_b,k)*k!)*args_list[4]^(1-sol_b+k), k, sol_b-sol_a, inf),
                f_2: factor * f_2,
                /*  we have -f_2 * diff(f_1, x) = (1-sol_b)! * (1-sol_b)! / pochhammer(sol_b-sol_a, 1-sol_b) * x^(-sol_b) + o(x^(-sol_b)) as x->0  */
                w0: (1-sol_b)! * (1-sol_b)! / pochhammer(sol_b-sol_a, 1-sol_b) * x^(-sol_b) * exp(x),
                wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
                )
            else
                (
                /*  https://dlmf.nist.gov/13.2.E30  */
                f_2: sum((1-sol_b)!*(k-1)!*args_list[4]^(1-sol_b-k)/((1-sol_b-k)!*pochhammer(sol_b-sol_a,k)), k, 1, 1-sol_b)
                    -sum(pochhammer(1+sol_a-sol_b,k)*args_list[4]^(1-sol_b+k)/(pochhammer(2-sol_b,k)*k!)
                        *(log(args_list[4])+psi[0](sol_a-sol_b+k+1)-psi[0](1+k)-psi[0](2-sol_b+k)), k, 0, inf),
                f_2: factor * f_2,
                /*  we have -f_2 * diff(f_1, x) = (1-sol_b)! * (1-sol_b)! / pochhammer(sol_b-sol_a, 1-sol_b) * x^(-sol_b) + o(x^(-sol_b)) as x->0  */
                w0: (1-sol_b)! * (1-sol_b)! / pochhammer(sol_b-sol_a, 1-sol_b) * x^(-sol_b) * exp(x),
                wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
                )
            )
        /* sol_b not integer  */
        else
            (
            /*  https://dlmf.nist.gov/13.2.E26  */
            f_1: factor * hypergeometric([sol_a], [sol_b], args_list[4]),
            f_2: factor * args_list[4]^(1-sol_b) * hypergeometric([ratsimp(1+sol_a-sol_b)], [ratsimp(2-sol_b)], args_list[4]),
            /*  https://dlmf.nist.gov/13.2.E33  */
            w0 : sin(%pi*sol_b)*exp(x)*gamma(sol_b)*gamma(2-sol_b)/(%pi*x^sol_b),
            wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )
        ),
    answer: write_solution(args_list, dep_var, x, f_1, f_2, wronski, answer),
    return(answer)
    )$


solve_bessel_or_kummer(args_list, h_function, dep_var, x, reg_sing, irreg_sing) := block(
    [answer, f_1, f_2, w0, wronski, x0, nm, det_data, coeffA, coeffB, coeffC, exponent, k, alpha, sol_a, sol_b, sol_b_2, phi, phi_inv, m_function_new, h_function_bessel, h_function_kummer, h_function_airy, args_new, factor],
    answer: [],
    if length(reg_sing)#1 or length(irreg_sing)#1 then 
        (
        print("Error: The function solve_bessel_or_kummer should only be called in case of 1 regular singular point and 1 irregular singular point.")
        )
    else if (rhs(reg_sing[1])=inf and rhs(irreg_sing[1])=inf) or (rhs(reg_sing[1])#inf and rhs(irreg_sing[1])#inf) then
        (
        print("Error: The function solve_bessel_or_kummer only works if infinity is either a regular or an irregular singular point.")
        )
    else
        (
        if rhs(reg_sing[1])#inf and rhs(irreg_sing[1])=inf then
            (
            x0: rhs(reg_sing[1])
            )
        else
            (
            x0: rhs(irreg_sing[1])
            ),
        nm: expand(ratsimp(args_list[3]*4*(x-x0)^2)),
        nm: expand(subst(x+x0,x,nm)),
        /* check whether nm(x) is of the form coeffA*(x-x0)^(2*exponent) + coeffB*(x-x0)^exponent + coeffC */
        det_data: detect_power_2(nm,x),
        coeffA: det_data[1],
        coeffB: det_data[2],
        coeffC: det_data[3],
        exponent: det_data[4],
        if diff(nm,x)=0 or exponent=inf or ratsimp(exponent)=0 then 
            (
            print("Error: The function solve_bessel_or_kummer only works if Q(x) is of the form P(x)/(4*(x-x0)^2) where P is of the form A*(x-x0)^(2*k)+B*(x-x0)^k+C with k not zero and A, B not both zero.")
            )
        else
            (
            if ratsimp(coeffB)=0 then
                /* in this case nm = coeffA * (x-x0)^exponent + coeffC   */
                (
                /*-------------------------------------------------------
                check for Airy function: q0(x)=(-4*alpha^3*k^2*(x-x0)^(3*k)+1-k^2)/(4*(x-x0)^2)
                and thus nm(x)=-4*alpha^3*k^2*(x-x0)^(3*k)+1-k^2
                --------------------------------------------------------*/
                k: exponent/3,
                if 1-k^2=coeffC then
                    (
                    alpha: (-coeffA/(4*k^2))^(1/3),
                    h_function_airy: 1,
                    phi: alpha*(x-x0)^k,
                    phi_inv: (x/alpha)^(1/k) + x0,
                    args_new: inv_transform(args_list, phi, phi_inv, x),
                    m_function_new: diff(args_new[4], x)^(1/2) * subst(args_new[4], x, h_function_airy),
                    print("The solutions can be expressed in terms of Airy functions."),
                    factor: h_function / m_function_new,
                    f_1: factor * airy_ai(args_new[4]),
                    f_2: factor * airy_bi(args_new[4]),
                    /* Wronskian for Airy functions see https://dlmf.nist.gov/9.2.E7  */
                    w0: 1 / %pi,
                    wronski: factor^2 * diff(args_new[4], x) * subst(args_new[4], x, w0),
                    answer: write_solution(args_new, dep_var, x, f_1, f_2, wronski, answer)
                    ),
                /*---------------------------------------------------------
                Bessel function B_(sol_a)(alpha*(x-x0)^k): nm(x) = 4*alpha^2*k^2*(x-x0)^(2*k) + 1-4*a^2*k^2
                and exponent=2*k
                ---------------------------------------------------------*/
                alpha: ratsimp(radcan(sqrt(ratsimp(coeffA)))/exponent),
                sol_a: ratsimp(radcan(sqrt(ratsimp(1-coeffC)))/exponent),
                h_function_bessel: x^(-1/2),
                phi: alpha*(x-x0)^(exponent/2),
                phi_inv: (x/alpha)^(2/exponent) + x0,
                args_new: inv_transform(args_list, phi, phi_inv, x),
                m_function_new: diff(args_new[4], x)^(1/2) * subst(args_new[4], x, h_function_bessel),
                answer: bessel_fund_system(args_new, h_function/m_function_new, dep_var, x, sol_a, answer),
                /* due to different notational conventions for exponent in the cases coeffB=0 and coeffB#0
                we need to divide exponent by 2 for the expression in terms of confluent hypergeometric functions*/
                exponent: exponent/2
                ),
            /*---------------------------------------------------------------
            expression in terms of confluent hypergeometric functions
            If coeffB#0 then nm = coeffA * x^(2*exponent) + coeffB * x^exponent + coeffC.
            Equations to solve for alpha, sol_a, sol_b:
            -alpha^2*exponent^2=coeffA,
            (2*sol_b-4*sol_a)*alpha*exponent^2=coeffB,
            1-(sol_b-1)^2*exponent^2=coeffC
            ---------------------------------------------------------------*/
            alpha: ratsimp(radcan(sqrt(ratsimp(-coeffA)))/exponent),
            /* sol_b, sol_b_2 two solutions to a quadratic equation */
            sol_b: ratsimp(1+radcan(sqrt(ratsimp(1-coeffC)))/exponent),
            sol_b_2: ratsimp(1-radcan(sqrt(ratsimp(1-coeffC)))/exponent),
            /* first consider sol_b */
            sol_a: ratsimp(sol_b/2-coeffB/(4*alpha*exponent^2)),
            h_function_kummer: x^(-sol_b/2)*exp(x/2),
            phi: alpha*(x-x0)^exponent,
            phi_inv: (x/alpha)^(1/exponent)+x0,
            args_new: inv_transform(args_list, phi, phi_inv, x),
            m_function_new: diff(args_new[4],x)^(1/2) * subst(args_new[4], x, h_function_kummer),
            answer: kummer_fund_system(args_new, h_function/m_function_new, dep_var, x, sol_a, sol_b, answer),
            /* consider also sol_b_2 */
            if ratsimp(sol_b_2-sol_b)#0 then
                (
                sol_a: ratsimp(sol_b_2/2-coeffB/(4*alpha*exponent^2)),
                h_function_kummer: x^(-sol_b_2/2)*exp(x/2),
                m_function_new: diff(args_new[4],x)^(1/2) * subst(args_new[4], x, h_function_kummer),
                answer: kummer_fund_system(args_new, h_function/m_function_new, dep_var, x, sol_a, sol_b_2, answer)
                )
            
            )
        ),
    return(answer)
    )$


solve_bessel_or_kummer_mobius_or_spheroid(args_list, h_function, dep_var, x, reg_sing, irreg_sing) := block(
    [beta, delta, phi, phi_inv, args_new, reg_sing_new, irreg_sing_new, answer],
    answer: [],
    if length(reg_sing)#2 or length(irreg_sing)#1 then 
        (
        print("Error: The function solve_bessel_or_kummer_mobius_or_spheroid should only be called in case of 2 regular singular points and 1 irregular singular point.")
        )
    else if reg_sing[1]=reg_sing[2] or reg_sing[1]=irreg_sing[1] or reg_sing[2]=irreg_sing[1] then 
        (
        print("Error: The function solve_bessel_or_kummer_mobius_or_spheroid only works if the equation has 3 distinct singular points.")
        )
    else if rhs(reg_sing[1])#inf and rhs(reg_sing[2])#inf and rhs(irreg_sing[1])#inf then
        (
        print("Error: The function solve_bessel_or_kummer_mobius_or_spheroid only works if infinity is a singular point.")
        )
    else if rhs(irreg_sing[1])=inf then
        (
        answer: solve_spheroid(args_list,h_function,dep_var,x,reg_sing)
        )
    /*-----------------------------------------------------------------------------------------------
    In the remaining case exactly one element in reg_sing is x=inf
    -----------------------------------------------------------------------------------------------*/
    else
        (
        if rhs(reg_sing[2])=inf then
            (
            beta: rhs(reg_sing[1])
            )
        else
            (
            beta: rhs(reg_sing[2])
            ),
        delta: rhs(irreg_sing[1]),
        phi: (x-beta)/(x-delta),
        phi_inv: (delta*x-beta)/(x-1),
        args_new: inv_transform(args_list, phi, phi_inv, x),
        reg_sing_new: [x=0],
        irreg_sing_new: [x=inf],
        answer: solve_bessel_or_kummer(args_new,h_function,dep_var,x,reg_sing_new,irreg_sing_new)
        ),
    return(answer)
    )$


solve_spheroid(args_list, h_function, dep_var, x, reg_sing) := block(
    [args_new, answer, f_1, f_2, wronski, dummy, phi, phi_inv, q0_mult, sol_q, sol_a, sol_b, sol_c, m_function_mathieu, m_function_new, h_function_mathieu, h_function_spheroid],
    answer: [],
    /*-----------------------------------------------------------------------------------------------
    We order the singular points increasingly in their real part
    -----------------------------------------------------------------------------------------------*/
    if imagpart(rhs(reg_sing[2])-rhs(reg_sing[1]))=0 and realpart(rhs(reg_sing[2]))>realpart(rhs(reg_sing[1])) then
        (
        dummy: reg_sing[2],
        reg_sing[2]: reg_sing[1],
        reg_sing[1]: dummy
        ),
    phi: (2*x-rhs(reg_sing[2])-rhs(reg_sing[1]))/(rhs(reg_sing[1])-rhs(reg_sing[2])),
    phi_inv: (rhs(reg_sing[1])-rhs(reg_sing[2]))*x/2+(rhs(reg_sing[1])+rhs(reg_sing[2]))/2,
    args_new: inv_transform(args_list, phi, phi_inv, x),
    q0_mult: expand(ratsimp(args_new[3]*4*(x-1)^2*(x+1)^2)),
    if diff(q0_mult,x,5)#0 then
        (
        print("The equation cannot be reduced to a spheroidal wave equation.")
        )
    else if coeff(q0_mult,x,3)#0 or coeff(q0_mult,x,1)#0 then
        (
        print("The equation cannot be reduced to a spheroidal wave equation.")
        )
    else
        (
        sol_q: coeff(q0_mult,x,4)/16,
        sol_a: ratsimp((coeff(q0_mult,x,0)-8*sol_q-2)/4),
        if asksign(ratsimp(coeff(q0_mult,x,2)+24*sol_q+4*sol_a-1))=zero then
            (
            print("Let mathieu_1(a,q,.), mathieu_2(a,q,.) be a fundamental system to Mathieu's equation"),
            print("f''(x)+(a-2*q*cos(2*x))*f(x)=0"),
            /*print("(1-x^2)*f''(x)-x*f'(x)+(a+2*q-4*q*x^2)*f(x)=0"),*/
            /*print("with a=", sol_a, "and q=", sol_q),*/
            h_function_mathieu: (1-x^2)^(-1/4),
            m_function_mathieu: diff(args_new[4],x)^(1/2)*subst(args_new[4],x,h_function_mathieu),
            print("The solutions can be given in terms of Mathieu functions."),
            f_1: h_function/m_function_mathieu*mathieu_1(sol_a,sol_q,acos(args_new[4])),
            f_2: h_function/m_function_mathieu*mathieu_2(sol_a,sol_q,acos(args_new[4])),
            wronski: tbc,
            answer: write_solution(args_new, dep_var, x, f_1, f_2, wronski, answer)
            ),
        sol_b: ratsimp(radcan(sqrt(4-coeff(q0_mult,x,4)-coeff(q0_mult,x,2)-coeff(q0_mult,x,0)))/2),
        sol_c: ratsimp(coeff(q0_mult,x,0)/4-sol_b-1),
        h_function_spheroid: (1-x^2)^(-(sol_b+1)/2),
        m_function_new: diff(args_new[4],x)^(1/2)*subst(args_new[4],x,h_function_spheroid),
        print("The solutions can be given in terms of spheroidal wave functions."),
        f_1: h_function/m_function_new*spherodialwave_a(sol_b,sol_c,sol_q,args_new[4]),
        f_2: h_function/m_function_new*spherodialwave_b(sol_b,sol_c,sol_q,args_new[4]),
        wronski: tbc,
        answer: write_solution(args_new, dep_var, x, f_1, f_2, wronski, answer)
        ),
    return(answer)
    )$




solve_legendre_or_hypergeo(args_list, h_function, dep_var, x, reg_sing, irreg_sing, reg_mult) := block(
    [dummy, phi, phi_inv, p0_new_mult, p1_new_mult, sol_nu, sol_mu, h_function_legendre, m_function_new, args_new, nm, sol_a, sol_b, sol_c, coeff2, coeffA, coeffB, coeffC, h_function_hypergeo, answer],
    answer: [],
    if length(reg_sing)#3 or length(reg_mult)#3 or length(irreg_sing)#0 then 
        (
        print("Error: The function solve_legendre_or_hypergeo only works if the equation has 3 regular singular points and 0 irregular singular points.")
        )
    else if reg_sing[1]=reg_sing[2] or reg_sing[1]=reg_sing[3] or reg_sing[2]=reg_sing[3] then 
        (
        print("Error: The function solve_legendre_or_hypergeo only works if the equation has 3 distinct regular singular points.")
        )
    else if rhs(reg_sing[1])#inf and rhs(reg_sing[2])#inf and rhs(reg_sing[3])#inf then 
        (
        print("Error: The function solve_legendre_or_hypergeo only works if the equation has a regular singular point at infinity.")
        )
    /*-----------------------------------------------------------------------------------------------
    In the remaining case exactly one element in reg_sing is x=inf
    -----------------------------------------------------------------------------------------------*/
    else
        (
        /*---------------------------------------------------------------------------------------------
        We put x=inf as the first element of reg_sing 
        ---------------------------------------------------------------------------------------------*/
        if rhs(reg_sing[1])#inf then
            (
            for i: 2 thru 3 do
            if rhs(reg_sing[i])=inf then
                (
                dummy: reg_sing[1],
                reg_sing[1]: reg_sing[i],
                reg_sing[i]: dummy,
                dummy: reg_mult[1],
                reg_mult[1]: reg_mult[i],
                reg_mult[i]: dummy
                )
            ),
        /*-----------------------------------------------------------------------------------------------
        We order the singular points increasingly in their real part
        -----------------------------------------------------------------------------------------------*/
        if imagpart(rhs(reg_sing[3])-rhs(reg_sing[2]))=0 and realpart(rhs(reg_sing[2]))>realpart(rhs(reg_sing[3])) then
            (
            dummy: reg_sing[2],
            reg_sing[2]: reg_sing[3],
            reg_sing[3]: dummy,
            dummy: reg_mult[2],
            reg_mult[2]: reg_mult[3],
            reg_mult[3]: dummy
            ),
        phi: (2*x-rhs(reg_sing[2])-rhs(reg_sing[3]))/(rhs(reg_sing[3])-rhs(reg_sing[2])),
        phi_inv: (rhs(reg_sing[3])-rhs(reg_sing[2]))*x/2+(rhs(reg_sing[3])+rhs(reg_sing[2]))/2,
        args_new: inv_transform(args_list, phi, phi_inv, x),
        p0_new_mult: expand(ratsimp(args_new[2]*(x^2-1)^2)),
        p1_new_mult: expand(ratsimp(args_new[1]*(x^2-1))),
        /*---------------------------------------------------------------------------------------------
        Case 1: the transformed equation is in fact a Legendre equation.
        The solutions are Legendre functions without additional factors.
        ---------------------------------------------------------------------------------------------*/
        if coeff(p0_new_mult,x,1)=0 and p1_new_mult=2*x then
            (
            /* print("The transformed equation is a Legendre equation"), */
            sol_nu: ratsimp((-1+radcan(sqrt(ratsimp(1-4*coeff(p0_new_mult,x,2)))))/2),
            sol_mu: sqrt(ratsimp(-coeff(p0_new_mult,x,2)-coeff(p0_new_mult,x,0))),
            /*if askinteger(sol_nu)=yes and askinteger(sol_mu)=yes then 
                (*/
            answer: legendre_fund_system(args_new, 1, dep_var, x, sol_nu, sol_mu, answer)
            /*    )*/
            )
        /*---------------------------------------------------------------------------------------------
        Case 2: the transformed equation is not a Legendre equation.
        The solutions are possibly Legendre functions with additional factors.
        ---------------------------------------------------------------------------------------------*/
        else
            (
            nm: expand(ratsimp(args_new[3]*4*(x^2-1)^2)),
            if diff(nm,x,3)#0 then
                (
                print("Error: The function Q passed to solve_legendre_or_hypergeo does not have infinity as a regular singular point.")
                )
            else
                (
                coeffA: coeff(nm,x,2),
                coeffB: coeff(nm,x,1),
                coeffC: coeff(nm,x,0),
                if coeffB#0 then 
                    (
                    print("The equation cannot be reduced to a Legendre equation.")
                    )
                else
                    (
                    sol_nu: ratsimp((-1+radcan(sqrt(ratsimp(1-coeffA))))/2),
                    sol_mu: radcan(sqrt(1-coeffA/4-coeffC/4)),
                    /*if askinteger(sol_nu)=yes and askinteger(sol_mu)=yes then
                        (*/
                    h_function_legendre: 1/sqrt(x^2-1),
                    m_function_new: diff(args_new[4],x)^(1/2)*subst(args_new[4],x,h_function_legendre),
                    answer: legendre_fund_system(args_new, h_function/m_function_new, dep_var, x, sol_nu, sol_mu, answer)
                        /*)*/
                    )
                )
            ),
        phi: (x-rhs(reg_sing[2]))/(rhs(reg_sing[3])-rhs(reg_sing[2])),
        phi_inv: (rhs(reg_sing[3])-rhs(reg_sing[2]))*x+rhs(reg_sing[2]),
        args_new: inv_transform(args_list, phi, phi_inv, x),
        p0_new_mult: expand(ratsimp(args_new[2]*x*(1-x))),
        p1_new_mult: expand(ratsimp(args_new[1]*x*(1-x))),
        /*---------------------------------------------------------------------------------------------
        Case 1: the transformed equation is in fact a hypergeometric equation.
        The solutions are hypergeometric functions without additional factors.
        ---------------------------------------------------------------------------------------------*/
        if diff(p0_new_mult,x)=0 and diff(p1_new_mult,x,2)=0 then
            (
            /*print("The transformed equation is a hypergeometric equation"),*/
            sol_c: ratsimp(coeff(p1_new_mult,x,0)),
            coeff2: coeff(p1_new_mult,x,1),
            sol_a: ratsimp((-1-coeff2-radcan(sqrt(expand((1+coeff2)^2+4*p0_new_mult))))/2),
            sol_b: ratsimp((-1-coeff2+radcan(sqrt(expand((1+coeff2)^2+4*p0_new_mult))))/2),
            answer: hypergeo_fund_system(args_new, 1, dep_var, x, sol_a, sol_b, sol_c, answer)
            )
        /*---------------------------------------------------------------------------------------------
        Case 2: the transformed equation is not a hypergeometric equation.
        The solutions are hypergeometric functions with additional factors.
        ---------------------------------------------------------------------------------------------*/
        else
            (
            nm: expand(ratsimp(args_new[3]*4*x^2*(x-1)^2)),
            if diff(nm,x,3)#0 then
                (
                print("Error: The function Q passed to solve_legendre_or_hypergeo does not have infinity as a regular singular point.")
                )
            else
                (
                coeffA: coeff(nm,x,2),
                coeffB: coeff(nm,x,1),
                coeffC: coeff(nm,x,0),
                /*---------------------------------------------------
                In some cases sign problems can arise in the following if the radicand contains variables.
                This is due to the behavior of the function radcan which I don't fully understand.
                -------------------------------------------------------*/
                sol_c: radcan(1+sqrt(ratsimp(1-coeffC))),
                sol_a: radcan((sol_c+sqrt(ratsimp(1-coeffA-coeffB-coeffC))+sqrt(ratsimp(1-coeffA)))/2),
                sol_b: radcan(sol_a-sqrt(ratsimp(1-coeffA))),
                if ratsimp(1-(sol_a-sol_b)^2-coeffA)#0 or ratsimp((1-sol_a-sol_b)*sol_c*2+4*sol_a*sol_b+coeffB)#0 or ratsimp(sol_c^2-2*sol_c+coeffC)#0 then
                    (
                    sol_a: radcan((sol_c-sqrt(ratsimp(1-coeffA-coeffB-coeffC))-sqrt(ratsimp(1-coeffA)))/2),
                    sol_b: radcan(sol_a-sqrt(ratsimp(1-coeffA))),
                    if ratsimp(1-(sol_a-sol_b)^2-coeffA)#0 or ratsimp((1-sol_a-sol_b)*sol_c*2+4*sol_a*sol_b+coeffB)#0 or ratsimp(sol_c^2-2*sol_c+coeffC)#0 then
                        (
                        print("Error: The function solve_hypergeo could not determine the coefficients")
                        )
                    else
                        (
                        /*print("Determination of a, b, c successful"),*/
                        h_function_hypergeo: exp(-(sol_c*log(x)+ratsimp(1+sol_a+sol_b-sol_c)*log(x-1))/2),
                        m_function_new: diff(args_new[4],x)^(1/2)*subst(args_new[4],x,h_function_hypergeo),
                        answer: hypergeo_fund_system(args_new, h_function/m_function_new, dep_var, x, sol_a, sol_b, sol_c, answer)
                        )
                    )
                else
                    (
                    /*print("Determination of a, b, c successful"),*/
                    h_function_hypergeo: exp(-(sol_c*log(x)+ratsimp(1+sol_a+sol_b-sol_c)*log(x-1))/2),
                    m_function_new: diff(args_new[4],x)^(1/2)*subst(args_new[4],x,h_function_hypergeo),
                    answer: hypergeo_fund_system(args_new, h_function/m_function_new, dep_var, x, sol_a, sol_b, sol_c, answer)
                    )
                )
            )
        ),
    return(answer)
    )$



/* see https://dlmf.nist.gov/14.2 */
legendre_fund_system(args_list, factor, dep_var, x, sol_nu, sol_mu, answer) := block(
    [f_1, f_2, w0, wronski],
    /* if realpart(nu)<-1/2 replace nu by -nu-1 */
    if asksign(realpart(sol_nu)+1/2)=neg then
        (
        sol_nu: -sol_nu-1
        ),
    if asksign(realpart(sol_mu))=neg then
        (
        sol_mu: -sol_mu
        ),
    /* thus we may assume nu + mu not in {-1, -2, ...}  */
    print("The solutions can be expressed in terms of associated Legendre functions"),
    print("of degree ", sol_nu, " and order ", sol_mu),
    print("(a fundamental system has not been implemented in all cases)."),
    if testmode#true and asksign(cabs(sol_nu)-cabs(sol_mu))#neg then
        (
        f_1: factor*assoc_legendre_p(sol_nu,sol_mu,args_list[4]),
        f_2: factor*assoc_legendre_q(sol_nu,sol_mu,args_list[4]),
        /* for w0 see Abramowitz and Stegun, Equation 8.1.8, page 333, without the factor exp(i*mu*pi) */
        /* it seems that assoc_legendre_p and assoc_legendre_q are implemented in Maxima without this factor */
        w0: 2^(2*sol_mu)*gamma((sol_nu+sol_mu+2)/2)*gamma((sol_nu+sol_mu+1)/2)/(gamma((sol_nu-sol_mu+2)/2)*gamma((sol_nu-sol_mu+1)/2)*(1-x^2)),
        wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0),
        answer: write_solution(args_list, dep_var, x, f_1, f_2, wronski, answer)
        ),
    /*print("(a fundamental system has not yet been implemented in all cases)"),
    if testmode#true and askinteger(sol_nu-sol_mu)=yes and asksign(sol_nu-sol_mu)=neg then
        (
        f_1: factor*assoc_legendre_p(sol_nu,sol_mu,args_list[4]),
        f_2: factor*assoc_legendre_p(sol_nu,sol_mu,-args_list[4]),
        wronski: tbc
        )
    else
        (
        f_1: factor*assoc_legendre_p(sol_nu,sol_mu,args_list[4]),
        f_2: factor*assoc_legendre_q(sol_nu,sol_mu,args_list[4]),
        wronski: tbc
        ),*/
    return(answer)
    )$



/*The function hypergeo_fund_system assigns
dep_var=%k1*factor*gauss_1(sol_a,sol_b,sol_c,args_list[4])+%k2*factor*gauss_2(sol_a,sol_b,sol_c,args_list[4])
where gauss_1 and gauss_2 are two linearly independent solutions of the hypergeometric equation.
The functions gauss_1 and gauss_2 are defined as in https://dlmf.nist.gov/15.10 */

hypergeo_fund_system(args_list, factor, dep_var, x, sol_a, sol_b, sol_c, answer) := block(
    [expand_hypergeometric: true, f_1, f_2, w0, wronski, f_list],
    print("The solutions can be expressed in terms of hypergeometric functions"),
    print("with parameters a =", sol_a, ", b =", sol_b, ", c =", sol_c),
    if testmode=true then
        (
        f_1: factor * gauss_a(sol_a,sol_b,sol_c,args_list[4]),
        f_2: factor * gauss_b(sol_a,sol_b,sol_c,args_list[4]),
        wronski: tbc
        )
    else if asksign(cabs(args_list[4])-1)=neg then
        (
        if sol_a=0 then
            (
            f_1: factor,
            f_2: factor * subst(args_list[4], s, integrate((1-s)^(sol_c-sol_b-1) * s^(-sol_c), s)),
            /* compute w0 by hand */
            w0: (1-x)^(sol_c-sol_b-1) * x^(-sol_c),
            wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )
        else if sol_b=0 then
            (
            f_1: factor,
            f_2: factor * subst(args_list[4], s, integrate((1-s)^(sol_c-sol_a-1) * s^(-sol_c), s)),
            /* compute w0 by hand */
            w0: (1-x)^(sol_c-sol_a-1) * x^(-sol_c),
            wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )
        else if askinteger(sol_c)=yes then
            (
            f_list: hypergeo_integer_case(dep_var, factor, sol_a, sol_b, sol_c, args_list[4], answer),
            f_1: f_list[1],
            f_2: f_list[2],
            wronski: f_list[3]
            )
        /*c not integer*/
        else
            (
            f_1: factor * hypergeometric([sol_a,sol_b], [sol_c], args_list[4]),
            f_2: factor * args_list[4]^(ratsimp(1-sol_c)) * hypergeometric([ratsimp(1+sol_a-sol_c), ratsimp(1+sol_b-sol_c)], [ratsimp(2-sol_c)], args_list[4]),
            /* https://dlmf.nist.gov/15.10.E3  */
            w0: (1-sol_c) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )
        )
    else if asksign(cabs(1-args_list[4])-1)=neg then
        (
        if sol_a=0 then
            (
            f_1: factor,
            f_2: factor * subst(args_list[4], s, integrate((1-s)^(sol_c-sol_b-1) * s^(-sol_c), s)),
            /* compute w0 by hand */
            w0: (1-x)^(sol_c-sol_b-1) * x^(-sol_c),
            wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )
        else if sol_b=0 then
            (
            f_1: factor,
            f_2: factor * subst(args_list[4], s, integrate((1-s)^(sol_c-sol_a-1) * s^(-sol_c), s)),
            /* compute w0 by hand */
            w0: (1-x)^(sol_c-sol_a-1) * x^(-sol_c),
            wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )
        else if askinteger(sol_c-sol_a-sol_b)=yes then
            (
            /* use Remark (d) on https://dlmf.nist.gov/15.10   */
            f_list: hypergeo_integer_case(dep_var, factor, sol_a, sol_b, 1+sol_a+sol_b-sol_c, 1-args_list[4], answer),
            f_1: f_list[1],
            f_2: f_list[2],
            wronski: f_list[3]
            )
        /*sol_c-sol_a-sol_b not integer*/
        else
            (
            f_1: factor * hypergeometric([sol_a,sol_b], [ratsimp(1+sol_a+sol_b-sol_c)], 1-args_list[4]),
            f_2: factor * (1-args_list[4])^(ratsimp(sol_c-sol_a-sol_b)) * hypergeometric([ratsimp(sol_c-sol_a), ratsimp(sol_c-sol_b)], [ratsimp(1+sol_c-sol_a-sol_b)],1-args_list[4]),
            /* https://dlmf.nist.gov/15.10.E5 */
            w0: (sol_a+sol_b-sol_c) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )
        )
    else if asksign(cabs(args_list[4])-1)=pos then
        (
        if sol_a=0 then
            (
            f_1: factor,
            f_2: factor * subst(args_list[4], s, integrate((1-s)^(sol_c-sol_b-1) * s^(-sol_c), s)),
            /* compute w0 by hand */
            w0: (1-x)^(sol_c-sol_b-1) * x^(-sol_c),
            wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )
        else if sol_b=0 then
            (
            f_1: factor,
            f_2: factor * subst(args_list[4], s, integrate((1-s)^(sol_c-sol_a-1) * s^(-sol_c), s)),
            /* compute w0 by hand */
            w0: (1-x)^(sol_c-sol_a-1) * x^(-sol_c),
            wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )
        else if askinteger(sol_a-sol_b)=yes then
            (
            /* use Remark (e) on https://dlmf.nist.gov/15.10 */
            f_list: hypergeo_integer_case(dep_var, factor*args_list[4]^(-sol_a), sol_a, sol_a-sol_c+1, sol_a-sol_b+1, 1/args_list[4], answer),
            f_1: f_list[1],
            f_2: f_list[2],
            wronski: f_list[3]
            )
        /*sol_a-sol_b not integer*/
        else
            (
            f_1: factor * args_list[4]^(-sol_a) * hypergeometric([sol_a, ratsimp(1+sol_a-sol_c)], [ratsimp(1+sol_a-sol_b)], 1/args_list[4]),
            f_2: factor * args_list[4]^(-sol_b) * hypergeometric([sol_b, ratsimp(1+sol_b-sol_c)], [ratsimp(1+sol_b-sol_a)], 1/args_list[4]),
            /* https://dlmf.nist.gov/15.10.E7  */
            w0: (sol_a-sol_b) * x^(-sol_c) * (x-1)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )
        )
    else
        (
        print("In this case a fundamental system for this equation in terms of hypergeometric functions has not yet been implemented."),
        f_1: factor*gauss_a(sol_a,sol_b,sol_c,args_list[4]),
        f_2: factor*gauss_b(sol_a,sol_b,sol_c,args_list[4]),
        wronski: tbc
        ),
    answer: write_solution(args_list, dep_var, x, f_1, f_2, wronski, answer),
    return(answer)
    )$


/* hypergeo_integer_case returns a fundamental system in the case that sol_c is an integer */
hypergeo_integer_case(dep_var, factor, sol_a, sol_b, sol_c, argument, answer) := block(
    [expand_hypergeometric: true, f_1, f_2, w0, wronski],
    
    /* sol_c in {1, 2, 3, ... } */
    
    if asksign(sol_c)=pos then
        (
        /*sol_a, sol_b in {0, -1, -2, ...} and sol_b >= sol_a:     https://dlmf.nist.gov/15.10.E10    */
        if askinteger(sol_a)=yes and asksign(sol_a)#pos and askinteger(sol_b)=yes and asksign(sol_b)#pos and asksign(sol_a-sol_b)#pos then
            (
            f_1: factor*hypergeometric([sol_a,sol_b],[sol_c],argument),
            f_2: factor*(hypergeometric([sol_a,sol_b],[sol_c],argument)*log(argument)
                -sum((sol_c-1)!*(k-1)!/((sol_c-k-1)!*pochhammer(1-sol_a,k)*pochhammer(1-sol_b,k))*(-argument)^(-k),k,1,sol_c-1)
                +sum(pochhammer(sol_a,k)*pochhammer(sol_b,k)*argument^k/(pochhammer(sol_c,k)*k!)
                *(psi[0](1-sol_a-k)+psi[0](1-sol_b-k)-psi[0](1+k)-psi[0](sol_c+k)),k,0,-sol_b)
                +(-1)^(-sol_b)*(-sol_b)!*sum((k-1+sol_b)!*pochhammer(sol_a,k)*argument^k/(pochhammer(sol_c,k)*k!), k, 1-sol_b, -sol_a)),
            /* w0 is a constant multiple of x^(1-sol_c)*(1-x)^(sol_c-sol_a-sol_b-1)  */
            /* The constant factor is determined by considering the term of order x^(-sol_c) in diff(f_2, x) */
            w0: (-1)^(sol_c-1)*(sol_c-1)!*(sol_c-1)!/(pochhammer(1-sol_a,sol_c-1)*pochhammer(1-sol_b,sol_c-1)) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(argument, x) * subst(argument, x, w0)
            )
        /*sol_a, sol_b in {0, -1, -2, ...} and sol_b < sol_a:     https://dlmf.nist.gov/15.10.E10   with sol_a, sol_b interchanged  */
        else if askinteger(sol_a)=yes and asksign(sol_a)#pos and askinteger(sol_b)=yes and asksign(sol_b)#pos and asksign(sol_a-sol_b)=pos then
            (
            f_1: factor*hypergeometric([sol_b,sol_a],[sol_c],argument),
            f_2: factor*(hypergeometric([sol_b,sol_a],[sol_c],argument)*log(argument)
                -sum((sol_c-1)!*(k-1)!/((sol_c-k-1)!*pochhammer(1-sol_b,k)*pochhammer(1-sol_a,k))*(-argument)^(-k),k,1,sol_c-1)
                +sum(pochhammer(sol_b,k)*pochhammer(sol_a,k)*argument^k/(pochhammer(sol_c,k)*k!)
                *(psi[0](1-sol_b-k)+psi[0](1-sol_a-k)-psi[0](1+k)-psi[0](sol_c+k)),k,0,-sol_a)
                +(-1)^(-sol_a)*(-sol_b)!*sum((k-1+sol_a)!*pochhammer(sol_b,k)*argument^k/(pochhammer(sol_c,k)*k!), k, 1-sol_a, -sol_b)),
            w0: (-1)^(sol_c-1)*(sol_c-1)!*(sol_c-1)!/(pochhammer(1-sol_a,sol_c-1)*pochhammer(1-sol_b,sol_c-1)) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(argument, x) * subst(argument, x, w0)
            )
        
        /*sol_a in {0, -1, -2, ...} and sol_b not in {sol_c-1, sol_c-2, ...., 0, -1, -2, ...}: https://dlmf.nist.gov/15.10.E9   */
        else if ( askinteger(sol_a)=yes and asksign(sol_a)#pos ) and ( askinteger(sol_b-sol_c)=no or asksign(sol_b-sol_c)#neg ) then
            (
            f_1: factor*hypergeometric([sol_a,sol_b],[sol_c],argument),
            f_2: factor*(hypergeometric([sol_a,sol_b],[sol_c],argument)*log(argument)
                -sum((sol_c-1)!*(k-1)!/((sol_c-k-1)!*pochhammer(1-sol_a,k)*pochhammer(1-sol_b,k))*(-argument)^(-k),k,1,sol_c-1)
                +sum(pochhammer(sol_a,k)*pochhammer(sol_b,k)*argument^k/(pochhammer(sol_c,k)*k!)
                *(psi[0](1-sol_a-k)+psi[0](sol_b+k)-psi[0](1+k)-psi[0](sol_c+k)),k,0,-sol_a)
                +(-1)^(-sol_a)*(-sol_a)!*sum((k-1+sol_a)!*pochhammer(sol_b,k)*argument^k/(pochhammer(sol_c,k)*k!), k, 1-sol_a, inf)),
            w0: (-1)^(sol_c-1)*(sol_c-1)!*(sol_c-1)!/(pochhammer(1-sol_a,sol_c-1)*pochhammer(1-sol_b,sol_c-1)) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(argument, x) * subst(argument, x, w0)
            )        
        /*sol_a not in {sol_c-1, sol_c-2, ...., 0, -1, ...} and sol_b in {0, -1, -2, ...}:   https://dlmf.nist.gov/15.10.E9  with sol_a, sol_b interchanged */
        else if ( askinteger(sol_b)=yes and asksign(sol_b)#pos ) and ( askinteger(sol_a-sol_c)=no or asksign(sol_a-sol_c)#neg ) then
            (
            f_1: factor*hypergeometric([sol_a,sol_b],[sol_c],argument),
            f_2: factor*(hypergeometric([sol_b,sol_a],[sol_c],argument)*log(argument)
                -sum((sol_c-1)!*(k-1)!/((sol_c-k-1)!*pochhammer(1-sol_b,k)*pochhammer(1-sol_a,k))*(-argument)^(-k),k,1,sol_c-1)
                +sum(pochhammer(sol_b,k)*pochhammer(sol_a,k)*argument^k/(pochhammer(sol_c,k)*k!)
                *(psi[0](1-sol_b-k)+psi[0](sol_a+k)-psi[0](1+k)-psi[0](sol_c+k)),k,0,-sol_b)
                +(-1)^(-sol_b)*(-sol_b)!*sum((k-1+sol_b)!*pochhammer(sol_a,k)*argument^k/(pochhammer(sol_c,k)*k!), k, 1-sol_b, inf)),
            w0: (-1)^(sol_c-1)*(sol_c-1)!*(sol_c-1)!/(pochhammer(1-sol_a,sol_c-1)*pochhammer(1-sol_b,sol_c-1)) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(argument, x) * subst(argument, x, w0)
            )
        
        /*sol_a, sol_b not in {sol_c-1, sol_c-2, ...., 0, -1, -2, ....}:   https://dlmf.nist.gov/15.10.E8    */
        else if ( askinteger(sol_a-sol_c)=no or asksign(sol_a-sol_c)#neg ) and ( askinteger(sol_b-sol_c)=no or asksign(sol_b-sol_c)#neg ) then
            (
            f_1: factor*hypergeometric([sol_a,sol_b],[sol_c],argument),
            f_2: factor*(hypergeometric([sol_a,sol_b],[sol_c],argument)*log(argument)
                -sum((sol_c-1)!*(k-1)!/((sol_c-k-1)!*pochhammer(1-sol_a,k)*pochhammer(1-sol_b,k))*(-argument)^(-k),k,1,sol_c-1)
                +sum(pochhammer(sol_a,k)*pochhammer(sol_b,k)*argument^k/(pochhammer(sol_c,k)*k!)
                *(psi[0](sol_a+k)+psi[0](sol_b+k)-psi[0](k+1)-psi[0](sol_c+k)), k, 0, inf)),
            w0: (-1)^(sol_c-1)*(sol_c-1)!*(sol_c-1)!/(pochhammer(1-sol_a,sol_c-1)*pochhammer(1-sol_b,sol_c-1)) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(argument, x) * subst(argument, x, w0)
            )
        
        /* sol_a in {1, ..., sol_c-1} or sol_b in {1, ..., sol_c-1}: see Remark (a) on https://dlmf.nist.gov/15.10 below Equation 15.10.7*/
        /* in particular, we have sol_c > 1 in this case  */
        else
            (
            f_1: factor*hypergeometric([sol_a,sol_b],[sol_c],argument),
            f_2: factor*argument^(ratsimp(1-sol_c))*hypergeometric([1+sol_a-sol_c,1+sol_b-sol_c],[2-sol_c],argument),
            /* https://dlmf.nist.gov/15.10.E3 */
            w0: (1-sol_c) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(argument, x) * subst(argument, x, w0)
            )
        )
    
    /* sol_c in {0, -1, -2, ...}*/
    /* According to Remark (c) on https://dlmf.nist.gov/15.10 multiply the previous solutions by z^(1-sol_c)           */
    /* and replace sol_a, sol_b, sol_c by 1+sol_a-sol_c, 1+sol_b-sol_c and 2-sol_c respectively */
    else
        (
        /*sol_a, sol_b in {sol_c-1, sol_c-2, sol_c-3, ...} and sol_b >= sol_a: https://dlmf.nist.gov/15.10.E10 with sol_a, sol_b, sol_c replaced */
        if askinteger(sol_a)=yes and asksign(sol_a-sol_c)=neg and askinteger(sol_b)=yes and asksign(sol_b-sol_c)=neg and asksign(sol_a-sol_b)#pos then
            (
            f_1: factor*argument^(1-sol_c)*hypergeometric([1+sol_a-sol_c,1+sol_b-sol_c],[2-sol_c],argument),
            f_2: factor*argument^(1-sol_c)*(hypergeometric([1+sol_a-sol_c,1+sol_b-sol_c],[2-sol_c],argument)*log(argument)
                -sum((1-sol_c)!*(k-1)!/((1-sol_c-k)!*pochhammer(sol_c-sol_a,k)*pochhammer(sol_c-sol_b,k))*(-argument)^(-k),k,1,1-sol_c)
                +sum(pochhammer(1+sol_a-sol_c,k)*pochhammer(1+sol_b-sol_c,k)*argument^k/(pochhammer(2-sol_c,k)*k!)
                *(psi[0](sol_c-sol_a-k)+psi[0](sol_c-sol_b-k)-psi[0](1+k)-psi[0](2-sol_c+k)),k,0,sol_c-sol_b-1)
                +(-1)^(sol_c-sol_b-1)*(sol_c-sol_b-1)!*sum((k+sol_b-sol_c)!*pochhammer(1+sol_a-sol_c,k)*argument^k/(pochhammer(2-sol_c,k)*k!), k, sol_c-sol_b, sol_c-sol_a-1)),
            w0: (-1)^(1-sol_c)*(1-sol_c)!*(1-sol_c)!/(pochhammer(sol_c-sol_a,1-sol_c)*pochhammer(sol_c-sol_b,1-sol_c)) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(argument, x) * subst(argument, x, w0)
            )
        /* sol_a, sol_b in {sol_c-1, sol_c-2, sol_c-3, ...} and sol_b < sol_a: https://dlmf.nist.gov/15.10.E10 with sol_a, sol_b, sol_c replaced */
        /* and sol_a, sol_b interchanged */
        else if askinteger(sol_a)=yes and asksign(sol_a-sol_c)=neg and askinteger(sol_b)=yes and asksign(sol_b-sol_c)=neg and asksign(sol_a-sol_b)=pos then
            (
            f_1: factor*argument^(1-sol_c)*hypergeometric([1+sol_b-sol_c,1+sol_a-sol_c],[2-sol_c],argument),
            f_2: factor*argument^(1-sol_c)*(hypergeometric([1+sol_b-sol_c,1+sol_a-sol_c],[2-sol_c],argument)*log(argument)
                -sum((1-sol_c)!*(k-1)!/((1-sol_c-k)!*pochhammer(sol_c-sol_b,k)*pochhammer(sol_c-sol_a,k))*(-argument)^(-k),k,1,1-sol_c)
                +sum(pochhammer(1+sol_b-sol_c,k)*pochhammer(1+sol_a-sol_c,k)*argument^k/(pochhammer(2-sol_c,k)*k!)
                *(psi[0](sol_c-sol_b-k)+psi[0](sol_c-sol_a-k)-psi[0](1+k)-psi[0](2-sol_c+k)),k,0,sol_c-sol_a-1)
                +(-1)^(sol_c-sol_a-1)*(sol_c-sol_a-1)!*sum((k+sol_a-sol_c)!*pochhammer(1+sol_b-sol_c,k)*argument^k/(pochhammer(2-sol_c,k)*k!), k, sol_c-sol_a, sol_c-sol_b-1)),
            w0: (-1)^(1-sol_c)*(1-sol_c)!*(1-sol_c)!/(pochhammer(sol_c-sol_a,1-sol_c)*pochhammer(sol_c-sol_b,1-sol_c)) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(argument, x) * subst(argument, x, w0)
            )
        
        /* sol_a in {sol_c-1, sol_c-2, ...} and sol_b not in {0, -1, -2, ....}: https://dlmf.nist.gov/15.10.E9 with sol_a, sol_b, sol_c replaced */
        else if ( askinteger(sol_a)=yes and asksign(sol_a-sol_c)=neg ) and ( askinteger(sol_b)=no or asksign(sol_b)=pos ) then
            (
            f_1: factor*argument^(1-sol_c)*hypergeometric([1+sol_a-sol_c,1+sol_b-sol_c],[2-sol_c],argument),
            f_2: factor*argument^(1-sol_c)*(hypergeometric([1+sol_a-sol_c,1+sol_b-sol_c],[2-sol_c],argument)*log(argument)
                -sum((1-sol_c)!*(k-1)!/((1-sol_c-k)!*pochhammer(sol_c-sol_a,k)*pochhammer(sol_c-sol_b,k))*(-argument)^(-k),k,1,1-sol_c)
                +sum(pochhammer(1+sol_a-sol_c,k)*pochhammer(1+sol_b-sol_c,k)*argument^k/(pochhammer(2-sol_c,k)*k!)
                *(psi[0](sol_c-sol_a-k)+psi[0](1+sol_b-sol_c+k)-psi[0](1+k)-psi[0](2-sol_c+k)),k,0,sol_c-sol_a-1)
                +(-1)^(sol_c-sol_a-1)*(sol_c-sol_a-1)!*sum((k+sol_a-sol_c)!*pochhammer(1+sol_b-sol_c,k)*argument^k/(pochhammer(2-sol_c,k)*k!), k, sol_c-sol_a, inf)),
            w0: (-1)^(1-sol_c)*(1-sol_c)!*(1-sol_c)!/(pochhammer(sol_c-sol_a,1-sol_c)*pochhammer(sol_c-sol_b,1-sol_c)) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(argument, x) * subst(argument, x, w0)
            )
        /* sol_a not in {0, -1, -2, ....} and sol_b in {sol_c-1, sol_c-2, ...}: https://dlmf.nist.gov/15.10.E9 with sol_a, sol_b, sol_c replaced */
        /* and sol_a, sol_b interchanged */
        else if ( askinteger(sol_b)=yes and asksign(sol_b-sol_c)=neg ) and ( askinteger(sol_a)=no or asksign(sol_a)=pos ) then
            (
            f_1: factor*argument^(1-sol_c)*hypergeometric([1+sol_b-sol_c,1+sol_a-sol_c],[2-sol_c],argument),
            f_2: factor*argument^(1-sol_c)*(hypergeometric([1+sol_b-sol_c,1+sol_a-sol_c],[2-sol_c],argument)*log(argument)
                -sum((1-sol_c)!*(k-1)!/((1-sol_c-k)!*pochhammer(sol_c-sol_b,k)*pochhammer(sol_c-sol_a,k))*(-argument)^(-k),k,1,1-sol_c)
                +sum(pochhammer(1+sol_b-sol_c,k)*pochhammer(1+sol_a-sol_c,k)*argument^k/(pochhammer(2-sol_c,k)*k!)
                *(psi[0](sol_c-sol_b-k)+psi[0](1+sol_a-sol_c+k)-psi[0](1+k)-psi[0](2-sol_c+k)),k,0,sol_c-sol_b-1)
                +(-1)^(sol_c-sol_b-1)*(sol_c-sol_b-1)!*sum((k+sol_b-sol_c)!*pochhammer(1+sol_a-sol_c,k)*argument^k/(pochhammer(2-sol_c,k)*k!), k, sol_c-sol_b, inf)),
            w0: (-1)^(1-sol_c)*(1-sol_c)!*(1-sol_c)!/(pochhammer(sol_c-sol_a,1-sol_c)*pochhammer(sol_c-sol_b,1-sol_c)) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(argument, x) * subst(argument, x, w0)
            )
        
        /* sol_a, sol_b not in {0, -1, -2, ....}: https://dlmf.nist.gov/15.10.E8 with sol_a, sol_b, sol_c replaced */
        else if ( askinteger(sol_a)=no or asksign(sol_a)=pos ) and ( askinteger(sol_b)=no or asksign(sol_b)=pos ) then
            (
            f_1: factor*argument^(1-sol_c)*hypergeometric([1+sol_a-sol_c,1+sol_b-sol_c],[2-sol_c],argument),
            f_2: factor*argument^(1-sol_c)*(hypergeometric([1+sol_a-sol_c,1+sol_b-sol_c],[2-sol_c],argument)*log(argument)
                -sum((1-sol_c)!*(k-1)!/((1-sol_c-k)!*pochhammer(sol_c-sol_a,k)*pochhammer(sol_c-sol_b,k))*(-argument)^(-k),k,1,1-sol_c)
                +sum(pochhammer(1+sol_a-sol_c,k)*pochhammer(1+sol_b-sol_c,k)*argument^k/(pochhammer(2-sol_c,k)*k!)
                *(psi[0](1+sol_a-sol_c+k)+psi[0](1+sol_b-sol_c+k)-psi[0](k+1)-psi[0](2-sol_c+k)), k, 0, inf)),
            w0: (-1)^(1-sol_c)*(1-sol_c)!*(1-sol_c)!/(pochhammer(sol_c-sol_a,1-sol_c)*pochhammer(sol_c-sol_b,1-sol_c)) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(argument, x) * subst(argument, x, w0)
            )
        
        /* sol_a, sol_b in {sol_c, sol_c+1, ..., 0} see Remark (a) on https://dlmf.nist.gov/15.10 below Equation 15.10.7*/
        else
            (
            f_1: factor*hypergeometric([sol_a,sol_b],[sol_c],argument),
            f_2: factor*argument^(1-sol_c)*hypergeometric([1+sol_a-sol_c,1+sol_b-sol_c],[2-sol_c],argument),
            /* https://dlmf.nist.gov/15.10.E3 */
            w0: (1-sol_c) * x^(-sol_c) * (1-x)^(sol_c-sol_a-sol_b-1),
            wronski: factor^2 * diff(argument, x) * subst(argument, x, w0)
            )
        ),
    return([f_1, f_2, wronski])
    )$


solve_4_reg_sing(args_list, h_function, dep_var, x, reg_sing, irreg_sing) := block(
    [dummy, q0_mult, delta, beta, alpha, phi, phi_inv, reg_sing_new, reg_mult_new, args_new, answer],
    answer: [],
    if length(reg_sing)#4 or length(irreg_sing)#0 then
        (
        print("Error: The function solve_4_reg_sing only works for 4 regular singular points and 0 irregular singular points.")
        )
    else if reg_sing[1]=reg_sing[2] or reg_sing[1]=reg_sing[3] or reg_sing[1]=reg_sing[4] or reg_sing[2]=reg_sing[3] or reg_sing[2]=reg_sing[4] or reg_sing[3]=reg_sing[4] then
        (
        print("Error: The function solve_4_reg_sing only works if the equation has 4 distinct regular singular points.")
        )
    else
        (
        /*---------------------------------------------------------------------------------------------
        We put x=inf as the first element of reg_sing 
        ---------------------------------------------------------------------------------------------*/
        if rhs(reg_sing[1])#inf then
            (
            for i: 2 thru 4 do
            if rhs(reg_sing[i])=inf then
                (
                dummy: reg_sing[1],
                reg_sing[1]: reg_sing[i],
                reg_sing[i]: dummy
                )
            ),
        q0_mult: ratsimp(args_list[3]*(x-rhs(reg_sing[2]))^2*(x-rhs(reg_sing[3]))^2*(x-rhs(reg_sing[4]))^2),
        if diff(q0_mult,x,3)#0 then
            (
            print("Error: The function solve_4_reg_sing only works if the numerator of Q(x) is a polynomial function of degree at most 2.")
            )
        else
            (
            delta: rhs(reg_sing[2]),
            beta: rhs(reg_sing[3]),
            alpha: (delta-rhs(reg_sing[4]))/(beta-rhs(reg_sing[4])),
            phi: alpha*(x-beta)/(x-delta),
            phi_inv: (delta*x-alpha*beta)/(x-alpha),
            args_new: inv_transform(args_list, phi, phi_inv, x),
            reg_sing_new: [x=inf,x=0,x=1],
            reg_mult_new: [unknown,2,2],
            answer: solve_legendre_or_hypergeo(args_new,h_function,dep_var,x,reg_sing_new,irreg_sing,reg_mult_new)
            )
        ),
    return(answer)
    )$
    

/*Let f=h*g where g satisfies g''+q0*g=0.
Let f_new(x)=f(phi(x)) and write f_new=h_new*g_new
where g_new satisfies g_new''+q0_new*g_new=0.
The function transform computes h_new and q0_new*/

transform(h_function, q0, phi, x) := block(
    [q0_new, h_function_new],
    q0_new: ratsimp(radcan(subst(phi,x,q0)*diff(phi,x)^2)+diff(phi,x,3)/(2*diff(phi,x))-3*diff(phi,x,2)^2/(4*diff(phi,x)^2)),
    h_function_new: ratsimp(radcan(subst(phi,x,h_function)*diff(phi,x)^(1/2))),
    return([h_function_new, q0_new])
    )$


/*We define args_list: [p1, p0, q0, argument, inhom].
Let f satisfy f''+p1*f'+p0*f=0 and write f=h*g
where g satisfies g''+q0*g=0.
Let phi and f_new such that f_new(phi(x))=f(x) and let p1_new, p0_new, q0_new 
such that f_new satisfies f_new''+p1_new*f_new'+p0_new*f_new=0
Write f_new(x)=h_new(x)*g_new(x) where g_new satisfies g_new''+q0_new*g_new=0.
The function inv_transform computes p1_new, p0_new and q0_new.*/

inv_transform(args_list, phi, phi_inv, x) := block(
    [p1_new, p0_new, q0_new, argument_new],
    q0_new: (args_list[3]+3*diff(phi,x,2)^2/(4*diff(phi,x)^2)-diff(phi,x,3)/(2*diff(phi,x)))/diff(phi,x)^2,
    q0_new: ratsimp(radcan(subst(phi_inv,x,q0_new))),
    /*print("Change of variables: the new Q is",factor(q0_new)),*/
    argument_new: factor(ratsimp(subst(args_list[4],x,phi))),
    p1_new: (args_list[1]+diff(phi,x,2)/diff(phi,x))/diff(phi,x),
    p1_new: ratsimp(radcan(subst(phi_inv,x,p1_new))),
    p0_new: args_list[2]/diff(phi,x)^2,
    p0_new: ratsimp(radcan(subst(phi_inv,x,p0_new))),
    /* We don't transform inhom. */
    return([p1_new, p0_new, q0_new, argument_new, args_list[5]])
    )$


determine_factorization(g,x):=block([k, i, j, coeff, sing_points, mult],
    /*This function checks whether g(x) is of the form coeff*(x-z1)^k1*...*(x-zn)^kn.
    If this is not true then the function returns [inf, [], []].
    If this is true then the function returns [coeff, [z1,...,zn], [k1,...,kn]].*/
    /*There are some problems with this function.
    Maybe it can be reformulated in a more efficient way using e.g.
    args(factor(g))
    This should give a list of all factors of g which can then be processed individually.
    Some of the factors may be of the form (x-x0)^k, some of the form (x^a+b)^k, etc.
    */
    if diff(g,x)=0 then
        (
        return([g, [], []])
        )
    else
        (
        sing_points: solve(factor(ratsimp(g))=0,x),
        mult: multiplicities,
        if length(sing_points)=0 then
            (
            print("Error: Cannot determine the zeros of the function", factor(ratsimp(g))),
            return([inf, sing_points, mult])
            )
        else
            (
            k: ratsimp(diff(g,x)/g),
            for i: 1 thru length(sing_points) do
                (
                k: ratsimp(radcan(k*(x-rhs(sing_points[i]))))
                ),
            if diff(k,x,length(sing_points))#0 then
                (
                print("Error: Cannot factor the function", factor(ratsimp(g))),
                return([inf, sing_points, mult])
                )
            else
                (
                coeff: g,
                for i: 1 thru length(sing_points) do
                    (
                    mult[i]: ratsimp(radcan(subst(rhs(sing_points[i]),x,k))),
                    for j: 1 thru i-1 do
                        (
                        mult[i]: ratsimp(radcan(mult[i]/(rhs(sing_points[i])-rhs(sing_points[j]))))
                        ),
                    for j: i+1 thru length(sing_points) do
                        (
                        mult[i]: ratsimp(radcan(mult[i]/(rhs(sing_points[i])-rhs(sing_points[j]))))
                        ),
                    coeff: ratsimp(radcan(coeff/(x-rhs(sing_points[i]))^(mult[i])))
                    ),
                if diff(coeff,x)#0 then
                    (
                    print("Error: Cannot factor the function", factor(ratsimp(g))),
                    return([inf, sing_points, mult])
                    )
                else
                    (
                    return([coeff, sing_points, mult])
                    )
                )
            )
        )
    )$


detect_power_1(g,x):=block([h, points, coeffA, exponent, x0],
    /*This function checks whether g(x) is of the form A*(x-x0)^exponent with constants A, exponent, x0.
    If this is so then it returns [A, exponent, x0].
    If this is not so then it returns [inf, inf, inf].*/
    if diff(g,x)=0 then
        (
        return([g, 0, 0])
        )
    else
        (
        h: ratsimp(diff(g,x)/g),
        if diff(ratsimp(diff(h,x)/h^2),x)=0 then
            (
            points: solve(1/h=0,x),
            if length(points)=0 or length(points)>1 then
                (
                return([inf, inf, inf])
                )
            else
                (
                x0: rhs(points[1]),
                exponent: ratsimp((x-x0)*h),
                coeffA: ratsimp(g/(x-x0)^exponent),
                if diff(exponent,x)#0 or diff(coeffA,x)#0 then
                    (
                    return([inf, inf, inf])
                    )
                else
                    (
                    return([coeffA, exponent, x0])
                    )
                )
            )
        else
            (
            return([inf, inf, inf])
            )
        )
    )$



detect_power_2(g,x):=block([h, coeffA, coeffB, coeffC, exponent],
    /*------------------------
    This function determines whether g(x) is of the form A*x^(2*k)+B*x^k+C or A*x^(2*k)+C or B*log(x)+C
    Ouput:
    If g(x)=A*x^(2*k)+B*x^k+C with k not 0 and B not 0 then the function returns [A, B, C, k].
    If g(x)=A*x^(2*k)+C with k not 0 then the function returns [A, 0, C, 2*k].
    If g is constant then the function returns [0, 0, g, 0].
    If g(x)=B*log(x)+C then the function returns [0, B, C, 0].
    If g(x) is not of any of these forms then it returns [inf, inf, inf, inf].
    ---------------------------*/
    if diff(g,x)=0 then
        (
        /*g is constant*/
        return([0, 0, g, 0])
        )
    else
        (
        h: ratsimp(diff(x*diff(g,x),x)/diff(g,x)),
        if diff(h,x)=0 then 
            (
            if h=0 then
                (
                /*g(x)=B*log(x)+C*/
                return([0, x*diff(g,x), radcan(g-coeffB*log(x)), 0])
                )
            else
                (
                /*g(x)=A*x^(2*k)+C*/
                exponent: h,
                coeffA: diff(g,x)*x^(1-h)/h,
                coeffB: 0,
                coeffC: ratsimp(g-coeffA*x^h),
                return([coeffA, coeffB, coeffC, exponent])
                )
            )
        else 
            (
            /*g is neither of the form A*x^(2*k)+C nor of the form B*log(x)+C*/
            exponent: ratsimp(diff(diff(x*diff(x*diff(g,x),x),x)/diff(g,x),x)/(3*diff(h,x))),
            if exponent=0 or diff(exponent,x)#0 then 
                (
                /*g is neither of the form A*x^(2k)+B*x^k+C*/
                return([inf, inf, inf, inf])
                )
            else 
                (
                coeffA: ratsimp(diff(x^(1-exponent)*diff(g,x),x)*x^(1-exponent)/(2*exponent^2)),
                if diff(coeffA,x)#0 then 
                    (
                    /*g is neither of the form A*x^(2k)+B*x^k+C*/
                    return([inf, inf, inf, inf])
                    )
                else
                    (
                    coeffB: ratsimp(diff(g-coeffA*x^(2*exponent),x)*x^(1-exponent)/exponent),
                    if diff(coeffB,x)#0 then 
                        (
                        /*g is neither of the form A*x^(2k)+B*x^k+C*/
                        return([inf, inf, inf, inf])
                        )
                    else
                        (
                        coeffC: ratsimp(g-coeffA*x^(2*exponent)-coeffB*x^exponent),
                        if diff(coeffC,x)#0 then 
                            (
                            /*g is neither of the form A*x^(2k)+B*x^k+C*/
                            return([inf, inf, inf, inf])
                            )
                        else
                            (
                            return([coeffA, coeffB, coeffC, exponent])
                            )
                        )
                    )
                )
            )
        )
    )$


detect_power_4(g,x):=block([coeffA, coeffC, exponent, x0, h1, h2, det_data, data_new],
    /*------------------------
    This function determines whether g(x) is of the form A*(x-x0)^(2*k-2)+C*(x-x0)^(k-2)
    with constants A, C, k#0, x0.
    Ouput:
    If g(x)=A*(x-x0)^(2*k-2)+C*(x-x0)^(k-2) with k#0 then the function returns [A, C, k, x0].
    If g(x) is not of any of these forms then it returns [inf, inf, inf, inf].
    ---------------------------*/
    det_data: detect_power_1(g,x),
    if det_data[2]#inf then
        (
        return([det_data[1],0,(det_data[2]+2)/2,det_data[3]])
        )
    else if diff(g,x,3)=0 then
        (
        coeffA: coeff(expand(g),x,2),
        coeffC: ratsimp(coeff(expand(g),x,0)-coeff(expand(g),x,1)^2/(4*coeff(expand(g),x,2))),
        exponent: 2,
        x0: ratsimp(-coeff(expand(g),x,1)/(2*coeff(expand(g),x,2))),
        return([coeffA, coeffC, exponent, x0])
        )
    else
        (
        h1: ratsimp(diff(g,x)/g),
        h2: ratsimp(diff(num(h1),x)),
        det_data: detect_power_1(h2,x),
        if det_data[2]#inf then
            (
            h1: ratsimp(subst(x+det_data[3],x,g)),
            data_new: detect_power_5(h1,x),
            if data_new[1]=inf then
                (
                return([inf, inf, inf, inf])
                )
            else
                (
                return([data_new[1], data_new[2], data_new[3], det_data[3]])
                )
            )
        else
            (
            return([inf, inf, inf, inf])
            )
        )
    )$


detect_power_5(g,x):=block([h1, h2, h3, exponent, coeffA, coeffC],
    /*------------------------
    This function determines whether g(x) is of the form A*x^(2*k-2)+C*x^(k-2)
    with constants A, C#0 and k#0
    Ouput:
    If g(x)=A*x^(2*k-2)+C*x^(k-2) with k not 0 then the function returns [A, C, k].
    If g(x) is not of any of these forms then it returns [inf, inf, inf].
    ---------------------------*/
    if g=0 then
        (
        return([inf, inf, inf])
        )
    else
        (
        h1: ratsimp(x*diff(g,x)/g),
        if diff(h1,x)=0 then
            (
            return([inf, inf, inf])
            )
        else
            (
            h2: ratsimp(x*diff(x*diff(g,x),x)/g),
            h3: ratsimp(diff(h2,x)/diff(h1,x)),
            if diff(h3,x)#0 then
                (
                return([inf, inf, inf])
                )
            else
                (
                exponent: ratsimp((h3+4)/3),
                if exponent=0 or diff(exponent,x)#0 then 
                    (
                    return([inf, inf, inf])
                    )
                else
                    (
                    coeffA: ratsimp(diff(g*x^(2-exponent),x)*x^(1-exponent)/exponent),
                    coeffC: ratsimp(g*x^(2-exponent)-coeffA*x^exponent),
                    if diff(coeffA,x)#0 or diff(coeffC,x)#0 then
                        (
                        return([inf, inf, inf])
                        )
                    else
                        (
                        return([coeffA, coeffC, exponent])
                        )
                    )
                )
            )
        )
    )$

kummer_fund_system_old(args_list, factor, dep_var, x, sol_a, sol_b, answer) := block(
    [expand_hypergeometric: true, f_1, f_2, wronski],
    if testmode=true then
        (
        print("The solutions can be expressed in terms of confluent hypergeometric functions."),
        print("Solution in testmode:"),
        f_1: factor * kummer_m(sol_a, sol_b, args_list[4]),
        f_2: factor * kummer_u(sol_a, sol_b, args_list[4]),
        wronski: tbc
        )
    else if askinteger(sol_a)=yes and asksign(sol_a)#pos and askinteger(sol_b)=no then
        (
        print("The solutions can be expressed in terms of generalized Laguerre polynomials"),
        print("of degree ", -sol_a, " with parameter ", ratsimp(sol_b-1), " and confluent hypergeometric functions."),
        if askinteger(sol_b)=no then 
            (
            f_1: factor * gen_laguerre(-sol_a, ratsimp(sol_b-1), args_list[4]),
            f_2: factor * args_list[4]^(ratsimp(1-sol_b)) * kummer_m(ratsimp(1+sol_a-sol_b), ratsimp(2-sol_b), args_list[4]),
            /* https://dlmf.nist.gov/13.2.E33 */
            w0: sin(%pi*sol_b) * x^(-sol_b) * exp(x) / %pi * gamma(sol_b) * gamma(2-sol_b),
            /* gen_laguerre(-a,b-1,x) is equal to pochhammer(b,-a)/(-a)! * kummer_m(a,b,x) by   https://dlmf.nist.gov/18.11.E2 */
            wronski: pochhammer(sol_b,-sol_a)/(-sol_a)! * factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )
        /* the function gen_laguerre(n, alpha, x) returns a division by 0 error e.g. for n=1, alpha=-1.
        Thus the following part is deactivated (this is the case askinteger(sol_a)=yes and asksign(sol_a)#pos and askinteger(sol_b)=yes ) */
        /*else
            (
            f_1: factor * gen_laguerre(-sol_a, ratsimp(sol_b-1), args_list[4]),
            f_2: factor * exp(args_list[4]) * kummer_u(ratsimp(sol_b-sol_a), sol_b, -args_list[4]),    */
            /*   https://dlmf.nist.gov/13.2.E38  */
        /*    w0: (-1)^(sol_a-sol_b) * x^(-sol_b) * exp(x),  */
            /* gen_laguerre(-a,b-1,x) is equal to (-1)^(-a)/(-a)! * kummer_u(a,b,x) by   https://dlmf.nist.gov/18.11.E2 */
        /*    wronski: (-1)^(-sol_a)/(-sol_a)! * factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )  */
        )
    else
        (
        print("The solutions can be expressed in terms of confluent hypergeometric functions."),        
        if is(sol_a=sol_b)=true then
            (
            f_1: factor * exp(args_list[4]),
            f_2: factor * kummer_u(sol_a, sol_b, args_list[4]),
            /* https://dlmf.nist.gov/13.2.E34  */
            w0: -x^(-sol_b) * exp(x),
            wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
            )
        else
            (
            if askinteger(sol_b)=yes and asksign(sol_b)#pos then 
                (
                /*  https://dlmf.nist.gov/13.2.E36  */
                f_1: factor * args_list[4]^(1-sol_b) * kummer_m(ratsimp(1+sol_a-sol_b), ratsimp(2-sol_b), args_list[4]),
                f_2: factor * kummer_u(sol_a, sol_b, args_list[4]),
                w0 : -x^(-sol_b)*exp(x)*gamma(2-sol_b)/gamma(sol_a-sol_b+1),
                wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
                )
            else if sol_a=1/2 and sol_b=3/2 then
                (
                f_1: factor * erf(%i*sqrt(args_list[4]))/sqrt(args_list[4]),
                f_2: factor * args_list[4]^(-1/2),
                /* https://dlmf.nist.gov/13.2.E34 */
                w0: -x^(-sol_b) * exp(x) * gamma(sol_a) / gamma(sol_b),
                /* erf(%i*sqrt(x))/sqrt(x) is equal to (2/sqrt(%pi)) * kummer_m(1/2, 3/2, x)  
                and x^(-1/2) is equal to kummer_u(1/2, 3/2, x)  */
                wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0) * 2 / sqrt(%pi)
                )
            else
                (
                f_1: factor * kummer_m(sol_a, sol_b, args_list[4]),
                f_2: factor * kummer_u(sol_a, sol_b, args_list[4]),
                /* https://dlmf.nist.gov/13.2.E34 */
                w0: -x^(-sol_b) * exp(x) * gamma(sol_a) / gamma(sol_b),
                wronski: factor^2 * diff(args_list[4], x) * subst(args_list[4], x, w0)
                )
            )
        ),
    answer: write_solution(args_list, dep_var, x, f_1, f_2, wronski, answer),
    return(answer)
    )$