@Preben Alsholm 

My approch should always fail if the solution is differentiable only in (0,1], being a polynomial approximation with a low degree.
Kitonum's example y2 is not a C^2 solution in [0,1]. I also wonder whether it is in (0,1]. I suspect that it is! An analytic example would be instructive (maybe for a similar problem).

@Preben Alsholm 

I have considered the ide as being defined in the closed interval [0,1](not in (0,1], or (0,1)), i.e. y''(0) exists and is finite. This does not imply that y'' is bounded near 0, but the boundedness is not needed.
We may complicate the problem and assume continuity in [0,1] and differentiability in (0,1], but is our intention to write a scientific paper?

Added. We may also drop the continuity at 0 and consider the integral as being improper, etc.

@Carl Love 

Unfortunately it works for a single unknown.

If y2(x) would be a solution <= exp(x) as it seems to be, then

z(x) = exp(x) - y2(x) satisfies the homogeneous equation

and z(0)=0, z(1)=exp(1)-1.5,  z >= 0.
But for x=0, the homogeneous equation implies int(z(t),t=0..1)=0, which is impossible.

The moral. A procedure to verify the solution is mandatory for this method.
Also, suppose that the solution is unique. Then how are the needed boundary conditions determined?

@cuongtd 

The norm is irrelevant.
Why don't you come with an example?
Here is a very simple one.

J:=<a,1,2; b,0,a; 1,a,3>:
V:=<-1,1,0>:
solve(convert(J.V, list));
                         {a = 1, b = 0}

Now, just try to find another suitable vector V.

You did not provide enough information or an example for a potential useful answer.
As I said, J.V=0 will probably not work unless V is chosen extremely careful.

@cuongtd 

@cuongtd 

If the rank does not matter, use det(J)=0.
Your second option will produce an incomatible system for most V.

@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.)