@J4James Could provide an example of first assigning to the constants and then using dsolve and yet obtaing imaginary results?

@J4James Do you get imaginary results if you use dsolve after assigning to the constants?

@samiyare In the worksheet I posted on January 27 I used different notation. I thought that worksheet better than my previous attempts.
So to be clear about this I need to get the whole worksheet you are using. Otherwise I shall be fumbling in the dark. Instead of pasting the whole code you could upload it using the thick green arrow in the MaplePrimes editor.

@Axel Vogt Let me try an explanation.

Any linear homogeneous ode with constant coefficients can be written as p(D)y = 0, where p is a polynomial function and D is the differentiation operator.
For all lambda and t we have
p(D)(exp(lambda*t) ) = p(lambda)*(exp(lambda*t)

Differentiating both sides with respect to lambda k times we get for all t and lambda

p(D)(t^k*exp(lambda*t) ) = kth derivative of p(lambda)*(exp(lambda*t) wrt. lambda.

Using the Leibniz rule on the right we get a sum of terms involving p(lambda) and derivatives of p(lambda) up to order k. The coefficients are linearly independent functions of t.

Thus if for some lambda y(t) = t^k*exp(lambda*t) is a solution to p(D)y = 0, then it follows by the linear independence that p and its derivatives up to order k are zero at lambda. Thus lambda is a root of order at least k+1 in p.

The converse also follows from that equation.

@kiyanoosh I have mentioned that you have a serious problem to be dealt with when the coefficient of the leading derivative becomes zero. This indeed must happen, since f ' (eta) -> 1 as eta -> infinity and f(0) = 0 are requirements. Thus it follows that f(eta) -> infinity as eta -> infinity making it certain that the coefficient of f ''' (eta), which is 1-pi*f(eta)^2, must go from being initially 1 to -infinity.
Your eq1 has the form
(1-pi*f(eta)^2)*f ''' (eta) + rest = 0,

where rest contains f, f ', f ''.

Rewriting you have

f ''' (eta) = -rest / ((1-pi*f(eta)^2)

If f ''' is to remain finite at the one or more points (eta1) where the denominator vanishes then rest must also vanish. This opens up the possibility of infinitely many solutions, one of which consists of the solution up to eta1 and after that continued as a solution to the 2nd order ode rest = 0.
This last one was the one I used in the second worksheet above.
To illustrate with simple examples the obvious problems take a look at these examples:
restart;
ode1:=(1-t)*diff(x(t),t)+2*x(t)=0;
dsolve({ode1,x(0)=1});
p0:=rhs(%);
p1:=piecewise(t<1,(t-1)^2,0);
odetest(x(t)=p1,ode1);
plot([p0,p1],t=0..2,thickness=3,smartview=false); #2 solutions
ode2:=surd(x(t),3)*diff(x(t),t)+3*x(t)=0; #Using surd(x(t),3) instead of x(t)^(1/3)
dsolve(ode2);
dsolve({ode2,x(0)=1});
q:=t0->piecewise(t<1,(1-t)^3,t<t0,0,(t0-t)^3);
plot([seq(q(t0),t0=1..2,0.1)],t=0..3,-1..1); #Infinitely many solutions

I don't think that that is a bug. For it to be a bug it should be claimed somewhere that is can do that kind of thing.
I have never seen such a claim.

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.