(For the next two commands see the installation page.)
libname := libname, FileTools:-JoinPath(["maple","lib","dcfun.mla"],base=homedir):
with(dcfun):
We choose a radix $b$ and a linear Mahler operator $L$.
b := 3;
$$3$$
L := x^3*(1-x^3+x^6)*(1-x^7-x^10) * M^2 - (1 - x^28 - x^31 - x^37 - x^40) * M + x^6*(1+x)*(1-x^21-x^30);
$$x^{3} \left(x^{6}-x^{3}+1\right) \left(-x^{10}-x^{7}+1\right) M^{2}-\left(-x^{40}-x^{37}-x^{31}-x^{28}+1\right) M +x^{6} \left(1+x \right) \left(-x^{30}-x^{21}+1\right)$$
It has order $r = 2$ and degree $d = 40$.
r, d := degree(L, M), degree(L, x);
$$2, 40$$
Let us search for the solutions $y$ of the homogeneous equation $Ly = 0$.
LMOpSolve(L, x, M, b, 'free' = 'K');
$$\frac{K_{1,1}}{\sqrt{x}}-K_{1,1} \sqrt{x}+K_{1,1} x^{\frac{3}{2}}-K_{1,1} x^{\frac{5}{2}}+K_{1,2} x^{3}+K_{1,1} x^{\frac{7}{2}}-K_{1,2} x^{4}-K_{1,1} x^{\frac{9}{2}}+K_{1,2} x^{5}+K_{1,1} x^{\frac{11}{2}}-2 K_{1,2} x^{6}+O\! \left(x^{\frac{13}{2}}\right)$$
By default, the procedure $\operatorname{LMOpSolve}$ search for solutions in the in the difference ring $${\mathfrak D} = \bigoplus_{\begin{subarray}{c}\lambda\in{\bar {\mathbb K}}\\ \lambda\neq 0\end{subarray}} \ln(1/x)^{\log_b(\lambda)} {\bar {\mathbb K}}((x^{1/*})), $$ where $\mathbb K$ is the base field and ${\bar {\mathbb K}}((x^{1/*}))$ is the field of Puiseux series over the algebraic closure ${\bar {\mathbb K}}$, equipped with the difference operator $M$.
In this example, we find a plan of Puiseux series with coefficients in the field of rational numbers or its algebraic closure, depending on how we interpret the constants $K_{1,1}$ and $K_{1,2}$. The first subscript of the constants $K$ indicates that $\lambda = 1$ is used.
We can search for ramified rational solutions too.
LMOpSolve(L, x, M, b, 'free' = 'K', 'output' = 'ramratpoly');
$$\frac{K_{1}}{\sqrt{x}\, \left(1+x \right)}$$
The result partially shows the origin of the operator.
The previous resolutions are based on the Newton polygon of the operator $L$ and on transformations.
Let us draw the Newton polygon of $L$.
LMOpNewtonPolygon(L, 0, x, M, b);
$$$$
The lower polygon, in blue color, corresponds to the search for series solutions, while the upper polygon corresponds to polynomial solutions. The edges of the lower polygon have slopes $-3$ and $1/2$, which means that solutions with valuation series $3$ and $-1/2$ are expected. More precisely, these are the only possible valuations, but there may be no solution with these valuations.
The search for Puiseux series solutions is reduced to the search for formal series solutions of an auxiliary operator.
L1 := LMOpChange(L, {x = t^2, y(x) = t^(-1)*z(t)}, x, M, b);
$$-t^{14} \left(t^{2}+1\right) \left(t^{60}+t^{42}-1\right)+\left(t^{80}+t^{74}+t^{62}+t^{56}-1\right) M -\left(t^{12}-t^{6}+1\right) \left(t^{20}+t^{14}-1\right) M^{2}$$
LMOpSolve(L1, t, M, b, 'free' = 'K');
$$K_{1,1}-K_{1,1} t^{2}+K_{1,1} t^{4}-K_{1,1} t^{6}+K_{1,2} t^{7}+\mathrm{O}\! \left(t^{8}\right)$$
LMOpSolve(L, x, M, b, 'free' = 'K');
$$\frac{K_{1,1}}{\sqrt{x}}-K_{1,1} \sqrt{x}+K_{1,1} x^{\frac{3}{2}}-K_{1,1} x^{\frac{5}{2}}+K_{1,2} x^{3}+K_{1,1} x^{\frac{7}{2}}-K_{1,2} x^{4}-K_{1,1} x^{\frac{9}{2}}+K_{1,2} x^{5}+K_{1,1} x^{\frac{11}{2}}-2 K_{1,2} x^{6}+O\! \left(x^{\frac{13}{2}}\right)$$
The previous computations deal with linear equations. We can also search for first-order right-hand factors $M - u$ of $L$. This amounts to solving a non-linear equation, which we call the Riccati equation associated with $L$, by analogy with the differential case.
To write this equation, we make the coefficients of the operator $L$ explicit.
for k from 0 to r do ell[k] := coeff(L, M, k) od;
$$x^{3} \left(x^{6}-x^{3}+1\right) \left(-x^{10}-x^{7}+1\right)$$
Riccati_equation := add(ell[k]*mul(u(x^(b^j)), j = 0..k-1) , k = 0..r);
$$x^{6} \left(1+x \right) \left(-x^{30}-x^{21}+1\right)+\left(x^{40}+x^{37}+x^{31}+x^{28}-1\right) u \! \left(x \right)+x^{3} \left(x^{6}-x^{3}+1\right) \left(-x^{10}-x^{7}+1\right) u \! \left(x \right) u \! \left(x^{3}\right)$$
We process. Note that we are using the operator $L$, not the previous writing.
LMOpSolve(L, x, M, b, 'free' = 'K', 'output' = 'riccati');
$$\left\{\frac{1}{x \left(x^{2}-x +1\right)}, x^{6} \left(1+x \right)\right\}$$
We obtain two solutions parametrized by the projective line ${\mathbb P}^1(\mathbb K)$, where $\mathbb K$ is the field extension of $\mathbb Q$ we want to consider.
The search for the solutions of the Riccati equation is equivalent to the search for what we call Mahler hypergeometric solutions.
LMOpSolve(L, x, M, b, 'free' = 'K', 'output' = 'mhypergeom');
$$\left\{K_{1} x^{3} \left(\overset{\infty}{\underset{\textit{\_k4} =0}{\textcolor{gray}{\prod}}}\! \frac{1}{1+x^{3^{\textit{\_k4}}}}\right), \frac{K_{2}}{\sqrt{x}\, \left(1+x \right)}\right\}$$