@Preben Alsholm OK, it's collocation then :-)
But usually in collocation methods the number of points equals the number of coefficients and a system is solved.
In my case I took nx >> n, that's why I think that choosing "special" points should not matter.
(I was not concerned about efficiency).

@Preben Alsholm 

With fsolve, the situation is clear.
The problem is with the minization approach.
With your example, this approach should find the exact solution (being in the basis). Here it's only NLP's fault of not being able to minimize properly.

I have tested also another example:

ide:=x^2*(diff(y(x), x, x))+(71/12)*x*(diff(y(x), x))-(7/4)*y(x)+int(y(t), t = 0 .. 1);
exactsol:=x^(1/3)+1;

which is not in the basis x^(k/2).

Added. By "no collocation" I meant not special chosen points.

@Preben Alsholm 

I tested you example withot collocation points and I got |y - yexact| < 2*10^(-6)

(This is normal because the exact solution is - unfortunately - in the basis!).

I wonder about the benefit of using collocation points instead of equidistant ones in this approach
(because the minimization is used).

@Preben Alsholm 

I just wanted to do something similar.
An analytic example would be nice.

Added. If someone is interested, here is y2 (without Chebychev polynomials and collocation)

Digits:=15:
ide:=x^2*(diff(y(x), x, x))+50*x*(diff(y(x), x))-35*y(x) - (1-exp(x+1))/(x+1)-(x^2+50*x-35)*exp(x)-int(exp(x*t)*y(t),t=0..1):
n:=10: nx:=16:
p:=x->add( a[k]*x^(k/2),k=0..n): 
F:=unapply( eval(ide, [y=p]),x):
d:=evalf(add(F(k/nx)^2,k=1..nx)):
s:=Optimization:-NLPSolve(d, {p(0)=1, p(1)=1.5}, iterationlimit=20000):
yy:=eval(p(x),s[2]);   # approx solution
yy=
   1.+.128100672246220*sqrt(x)-4.96593663847107*x+21.9900671494260*x^(3/2)-
   61.2588096618634*x^2+119.676248631378*x^(5/2)-158.798014306812*x^3+
   142.246959241914*x^(7/2)-81.9327391785196*x^4+27.4658220286035*x^(9/2)-4.05169793790087*x^5
The plot is the same as Preben's.

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