Nous définissons dans un premier temps l'expression Q comme suit:
> Q := x^2 * diff(f(x),x,x) + x*diff(f(x),x) - p^2*f(x);
La substitution de f(x) par le monôme 
se fait par:
> subs(f(x)=x^k,Q);
Il ne reste plus qu'à évaluer cette dernière expression et à trouver les valeurs de k qui l'annulent.
> expand(");
> solve(",k);
Ce qui nous donne k=p ou k=-p. L'équation différentielle ordinaire
étant du second ordre, elle admet au plus deux degrés de liberté. Les
deux solutions 
et 
sont indépendantes, donc elles
forment une base pour l'ensemble des solutions. Donc la forme générale
des solutions de l'équation différentielle est une combinaison
linéaire des deux solutions indépendantes, i.e. 
. On
retrouve évidemment ce résultat avec l'instruction dsolve(Q,f(x)).
On pose:
> F : = A*x^p + B/x^p;
On introduit les conditions aux limites, puis on résoud le système en A et B, et enfin on réintroduit les solutions de A et B dans l'expression de F:
> solve({subs(x=R0,F)=1, subs(x=R0,diff(F,x))=0},{A,B});
> subs(",F);
ce qui aboutit au résultat final:
On aurait abouti au même résultat en utilisant la syntaxe suivante:
> dsolve({Q,f(R0)=1,D(f)(R0)=0},f(x));