@J4James With z=P/L^2*exp(-L*x), T(x) = U(z)  and bcs:=T(0)=1,T(infinity)=0;

we have that x = 0 corresponds to z = P/L^2.
But x = infinity corresponds to z = 0 only if L > 0.
If L < 0 then x = infinity corresponds to z = infinity, since then exp(-L*x) -> infinity as x -> infinity.
(I'm assuming P > 0 in this the latter case too. Otherwise x = infinity corresponds to -infinity).

@J4James I think you miscalculated some.
But the following does come up with a relatively compact answer. The boundary conditions for U are assuming that L > 0.

restart;
ODE:=(diff(T(x), x, x))+P*(S+a*(1-exp(-L*x))/L)*(diff(T(x), x))=0;
bcs:=T(0)=1,T(infinit)=0;
var:=z=P/L^2*exp(-L*x);
var2:=x=solve(var,x);
PDEtools:-dchange({var2,T(x)=U(z)},ODE,[z,U(z)]);
collect(%,diff,factor);
ode1:=isolate(%,diff(U(z),z,z));
dsolve({ode1,U(P/L^2)=1,U(0)=0});

It seems to me from the image you presented that the first x is right next to the first parenthesis (.
This means that there is a missing multiplication sign. The same appears to be true for the first y.

@J4James Yes indeed. The exponential exp(4.77*z) seems to be the underlying problem.

Try with infinity = 1 just to see what is going on.
I noticed that your L is in fact negative, so I also tried replacing L with -L1 in in analytical computation of the solution. The solutions returned in the two cases used different special functions. THe one with Whittaker performed better at the end.

restart:
ODE:=diff(T(z),z$2)+A1*(S-1/L+1/L*exp(-L*z))*diff(T(z),z)+A2*T(z)=0;
bcs:=T(0)=1,T(inf)=0;
ODE2:=subs(L=-L1,ODE);
sol2:=dsolve({ODE2,bcs});
sol1:=dsolve({ODE,bcs});
################
L:=(S-sqrt(S^2+4*alpha*epsilon))/(2*epsilon);
A1:=3*R*P/(epsilon2*(3*R+4));
A2:=3*n*R*P/(epsilon2*(3*R+4));
L1:=-L;
plot(eval(rhs(sol2),{epsilon=.01,epsilon2=.1,alpha=5,R=1,n=1,P=6,S=1,inf=1}),z=0..1);
plot(eval(rhs(sol1),{epsilon=.01,epsilon2=.1,alpha=5,R=1,n=1,P=6,S=1,inf=1}),z=0..1);
########Numerically
Digits:=15:
ODE;
BVP:=eval({ODE,bcs},{epsilon=.01,epsilon2=.1,alpha=5,R=1,n=1,P=6,S=1,inf=1});
sol:=dsolve(BVP,numeric,maxmesh=1024);
plots:-odeplot(sol,[z,T(z)],0..1);

@Luka You got infinitely many solutions with 4 free parameters, although one of them t4 must be positive as we noticed. With res being defined as above you can find the free parameters as the ones appearing in res as ti = ti:

lhs~(select(evalb,res));

I get {t3, t4, t6, t9}.

@J4James Well, the main problem is the unevaluated limit in the result. That is clearly due to the boundary condition at infinity. Replace it (as in the original question) by a symbol, say inf, then

bcs:=T(0)=1,T(inf)=0;

sol:=dsolve({ODE,bcs}):
plot(eval(rhs(sol),{epsilon=.01,epsilon2=.1,inf=5}),eta=0..5);

Whether or not the result contains imaginary numbers is irrelevant. Think of
I*(exp(I*t)-exp(-I*t)) for t real.

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