@Preben Alsholm 

I have corrected the approx solution, thank you for detecting the error.

An interesting fact is that the wrong solution was very close to the exact one.
Of course in the present version it is better.

@Carl Love 

Unfortunately I don't know. The existence/uniqueness of the solution should be studied theoretically, the numerical solutions of such problems are not easy.
I just tried to find a simple function with small residuals. A polynomial is convenient because the associated minimization problem is quadratic.

@Kitonum 

It is not a bug, it is the (generalized) sum of a divergent series (well known in mathematics) and it is documented, see ?evalf,Sum. It can be turned off using _EnvFormal:=false (or the formal option in Maple 2016).

@tomleslie 

should be return 1/(1-w)

@Thomas Richard 

- parameters seems to be documented only for numeric ODEs. What is its role here?

- Why without parameters one obtains:
   Error, (in dsolve) numeric exception: division by zero

- Without initial conditions, i.e.
dsolve(sys, {x(t),y(t)});
==> 4 solutions,  but the first two are not!

@Kitonum 

Nice solution, thumb up!
A bit simplified version:

sys:={diff(x(t),t) = -k/m*x(t)*sqrt(x(t)^2+y(t)^2),
           diff(y(t),t) = -k/m*y(t)*sqrt(x(t)^2+y(t)^2)}:
Polar:={x(t)=r(t)*cos(phi(t)), y(t)=r(t)*sin(phi(t))}:
syspolar:=simplify(eval(sys,Polar)) assuming r(t)>0:
solpolar:=dsolve(syspolar union {r(0)=r__0,phi(0)=phi__0},{r(t),phi(t)}):
sol:=eval(Polar, solpolar);

@ecterrab 

The answer is for the "non-strange" PDE where u(0,t) is replaced by u(x,t).

@Carl Love 

I was at first thinking at

A1:=map(eval, A, consts):
A2:=map(eval, A, c=66):

(which should work in Maple <10), and I was happy to see that map is not necessary in Maple 2016. (I did not see that OP uses Maple 15.)

@Christopher2222 

You can simply re-define the ranking of the loser(s) to 0.

@asa12 

I simply don't understand your maths and don't know about what solutiion you are talking.

The problem is very simple. GL(2,3) is (isomprphic to) a subgroup of S_8 generated by two permutations (those given by Generators(GL23)). So, it is also isomorphic to a group of permutation matrices. That's all.

@asa12 

Of course you can convert to matrices, but apriori there is none!

GL23 := GeneralLinearGroup(2, 3):
g:=Generators(GL23):
J:=LinearAlgebra:-IdentityMatrix(8):
J[convert(g[1],permlist,8)];  # this should be your aa;  

What solution? Having the generators, any element of the group is a product of them.
The sum of two matrices in GL is not generally in GL.

@asa12 

You were told that matrix representation is not implemented.
Actually, for GL(n,q) such a representation is obvious by definition and you may generate easily the matrices (al least when q is prime).

Representations as permutation groups are available only for GeneralLinearGroup and GeneralOrthogonalGroup, but for small parameters only (see the documentation).

@systemcode 

I do not know the algorithm in detail, the code is not simple, see:
showstat(fsolve);

@gkokovidis 

No, I have changed a bit the example such that Digits=400 is not enough.

@gkokovidis 

fsolve({(2*x+y+1)*exp(-(x-y-2)^4),   (3*x+2*y-1)*exp(-(x-y-1)^4)});