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)

@gawati2611

You must download (from Maple Application Center) and install the package. 

I think you are supposed to do the computations by hand (it's very easy) and then check them with Maple.
If you don't, you will never know linear algebra!

dharr 

@ReactionUra 

Some variables are mathematically negative; they are small in absolute value, so probably they can be considered 0.

Note that dharr's solution is similar; actually if you start his worksheet with Digits:=25, fsolve fails just because of those variables.

Explore(plot([[[0,0],[1,0],[1/3,2],[0,0]],[(12*s*(t-1)*(s-1)*xi^2+t)/(3+12*(t-1)*(s-1)*xi^2+12*(t-1)*(s-1)*xi),2*t/(1+4*(t-1)*(s-1)*xi^2+4*(t-1)*(s-1)*xi), xi=-infinity..infinity]]), parameters=[s=0. .. 1, t=0. .. 1]);