I tried looking at an initial value problem for this system, i.e replacing the conditions at t = 50 with conditions at t = 0,
x[1](0)=1.7*10^8, x[2](0)=0,x[3](0)=400,psi[1](0)=p1,psi[2](0)=p2,psi[3](0)=p3,
and used the parameters option in dsolve.
The behavior I saw while varying parameters left me believing that there is no chance of solving the bvp over an interval that long.
I succeeded only in solving the bvp on t = 0 .. 0.45, i.e. on an interval less than 100 times smaller.

@samiyare The shorter version I mentioned I shall test a little more, but it worked fine with N[bt]=0.4.
If it appears that it covers a greater range of parameter values without too much advance experimentation then I shall upload that version.
However, I think that these types of problems are tricky.

And an animation could be made:

eq1 := diff(y(x), x) = -(2*x*cos(y(x))+3*x^2*y(x))/(x^3-x^2*sin(y(x))-y(x));
s:=dsolve(eq1);
plots:-animate(plots:-implicitplot,[s, x=-4..4, y=-4..4, gridrefine=2],_C1=-5..5 ,trace=24);

@Carl Love The Maple choice of parentheses () for grouping as well as for use in function application has its drawbacks, but it does follow ordinary mathematical practice. The 2D interface has probably just added to the problems. Mathematica made another choice.

I'm not so certain that Maple should issue an error message. It has always accepted 'weird' input like sin(x)(x) or x(x).
It plots both and integrates both numerically since real or complex numbers are accepted as constant functions. Replacing x with 2 above as in sin(2)(2) or 2(2) gives us sin(2) and 2, respectively.

Belaboring some:
u:=(exp(x)-1)(1-exp(-x));
ustar:=(exp(x)-1)*(1-exp(-x));
g:=(exp(x)-1);
plot([u,g,ustar],x=0..1,color=[red,blue,green]);
eval~([u,g,ustar],x=ln(2));
map(int,[u,g,ustar],x=0..1);
evalf(%);
int(x(x),x=0..1);
evalf(Int(x(x),x=0..1));

@acer I didn't catch that one!

@ 
The warning you got usually means that the optimization method erroneously computed the gradient to be zero.
However, to test the problem I made up my own data from the known values of alpha = 0.0035, theta = 60.
Obviously this should make for a good solution, but my point was to see if it iterated or not.
It did and it did find the solution alpha = 0.0035, theta = 60.

restart;
with(Statistics):
rho[0]:=1.33;
T0:=6;
f:=proc(T,alpha,theta) evalf(rho[0]+alpha*(T-T0)*((T-T0)/theta)^5*(Int(x^5/(exp(x)-1)(1-exp(-x)) ,x = 0 .. theta/(T-T0)))) end proc;
temp:=<seq(7..25)>;
fT:=rcurry(f,.0035,60);
resistance:=fT~(temp);
NonlinearFit(f,temp,resistance,output=[parametervalues, residuals],parameterranges=[0.001..0.01,10..100],initialvalues=[0.0027,50]);

Note:
I also tried to perturb the data slightly. NonlinearFit came up with a result here too:
X := RandomVariable(Normal(0, 1e-4));
A := Sample(X, 19):
resistance2:=resistance+convert(A,Vector[column]);
NonlinearFit(f,temp,resistance2,output=[parametervalues, residuals],parameterranges=[0.001..0.01,10..100],initialvalues=[0.0027,50]);

It should be NonlinearFit (not NonLinearFit).
You should use := for assignments.
Do not use T[0] when T is in use on its own. Use T0 (or whatever) instead.
The model has an unevaluated integral. Maple cannot do the integral, not even when the upper limit is 1.
You may want to try the operator form using something like

f:=proc(T,alpha,theta)
  evalf(rho[0]+alpha*(T-T0)*((T-T0)/theta)^5*(Int(x^5/(exp(x)-1)(1-exp(-x)) ,x = 0 .. theta/(T-T0))))
end proc;
#Test of f
f(7,.0027,50);

We don't get much further without the data.

@lcrpn 
Here is a slightly shorter version:

expand(f(x));
collect(%,sin(x));
convert(%,phaseamp,x);
res:=simplify(%);
#Incidentally the argument to cos can be simplified:
w:=op([2,1,2],res);
#Getting an idea:
identify(evalf(w));
#Showing that w+3*Pi/20 = 0:
u:=w+3*Pi/20;
expand(tan(u)); #we know that -Pi/2 < w <0 so that -Pi/2 < u < Pi/2
simplify(%);

I don't know what is going on, but noticed that if a is given a concrete value from the start (instead of assuming on it) e.g. a:=1/3 or a:=1/2 there is no difference.

There are a couple of syntax problems:

A missing od (end do) after the line
W[i]:=unapply(%,x):                 #Orthonormality process

f[j](x) should be replaced by f[j] in the line
c[io][j]:=int(f[j](x)*W[io](x),x=0..1);

I strongly suggest replacing linalg[matrix] by Matrix. In fact for the code you have I would use Vector(2^J+2) instead of Matrix(2^J+2,1). However, I have only looked at the syntax of the problem as presented in the code.

To see output from the two double loops at the end set printlevel at 2:
printlevel:=2;

There is no file attached. Use the thick green arrow in the icon panel of the editor.

@mahmood180 By applying the laplace transform to the ode you get an ode for the laplace transform. I Using dsolve on that ode you get an integral that seems rather undoable.

So I only see iteration as a way to solve your new ode. That method would also solve your previous one.
What is needed is knowledge of x(t) at some interval of length a = 1/4 then you can work your way up.
This will be an exact solution.
Similar to the previous problem I'm here assuming that x(t) = 0 for t < 0 and x(0) = x0. This implies that x(t) = x0 on t = 0..a, thus we have a start.

I posted the iterative solution in a general form as an answer and have removed the remainder of this comment.

@ateaktree The upload of the file didn't succeed. I got an error 404.

@mahmood180 Certainly. If I think I can contribute in any way I shall be happy to.
MaplePrimes is the perfect forum.

@mahmood180 
There is obviously something I am not getting. It doesn't seem to me that your T1 is a solution.

restart;
ode:=diff(x(t),t)=x(t-1)+2*t;
#If x(t)=0 for t<0 and x(0)=x0 does it not follow that x(t)=t^2+x0 for 0 <= t < 1?
eval(ode,{x(t)=t^2+x0,x(t-1)=0}); #Certainly is a solution