@mmcdara 

Why do you think that there is "the" exact solution. Maybe the ode has more than one (standard uniqueness theorems do not apply here).

@nrussoniello 

Each symbol is replaced with a random number, e.g. s1=77, s2=-665 ...

Suppose that the entries of a 2x2 matrix M(x,y) are polynomials in x,y.
The rank of M(x,y) is <=1 iff det(M(x,y)) is identical 0. E.g. if det(M(x,y)) = x^2 - y + x*y^7 then rank(M(x,y)) = 2.

If follows that if det(M(78,55)) <>0 then rank(M(x,y))=2. 
Your case is similar.

The sitation would have been more complicated if the rank for the matrix with symbols replaced were not maximal.Then we must repeat the procedure and take the maximal ranks, but accepting now that the true rank could be higher (but with a small probability).

If you don't want to read at least a tutorial, the shortest path would be to post here a simplified version of your problem.
Probably someone will upload a worksheet with the solution, and you will be able to adapt it for similar tasks.

@Carl Love 

This order is not total (linear) and has nothing to do with "lexicographic". Where is it used? How is x[2] computed using this order?

x[n] seems to be the least positive interger not in { x[1],...,x[n-1] } which satisfies ...

Do you mean the lexicographic order in N? But then 123 < 2 < 3, so why is x[2] = 3?

It is not only old code, but you have also replaced the old keyword array with the new one Array in a few places. Revert these changes.

@mmcdara The first solution is very elegant; of course I like it. But it was probably provided by the author/committee or maybe by a very smart student (who deserves congratulations & a special prize, if he/she exists).

I don't know what experience you have in this field, but I was involved in these competitions
as a direct participant (a long time ago) and then as a team leader for my university students
(I have even proposed problems for the IMC contest).

Now I don't have such preoccupations but I still enjoy to occasionally watch these competitions
and find solutions (with or without Maple).
 

I think it's better to simplify first by hand.
Note first that expr = 0 for r>n (both terms are 0).
Then express everything w.r.t.
y = q^r;
x = q^n;
Y = product(y - q^i, i = 1 .. r - 1);
X = product(x - q^i, i = 1 .. r - 1);
and finally, simplify (or factor) the result (==> rational function in q,x,y,X,Y).

Maple simplifies telescoping sums:

@nm You are confusing the inverse of a function with multiplicative inverse.

@nm 

Of course I would not classify y = g(y') as d'Alembert, but it matches the pattern y = x*f(y') + g(y') for f = 0
and the solution is valid for f=0 too. Anyway, the classification is mainly didactic and Maple does non necessarily use the classical methods. So, your problem is more or less artificial.

@acer Ok, so, a bug in IntegrationTools:-Expand (2015).

@mthkvv 

modp1 is builtin, so its code is not accessible. You will need to be employed by Maplesoft :-)

@gawati2611 

https://www.maplesoft.com/applications/view.aspx?SID=33406

(it is the first result with any search engine)