From previous results the following guess for a solution should roughly have the correct shape:
GUESS:=[diff(u(eta), eta) = .12-.25*eta, T(eta) = eta*(0.7*eta-1), u(eta) = .12*eta-.125*eta^2, diff(T(eta), eta) = -1+1.4*eta, phi(eta) = .1*eta*(1.4-eta), phi0 = 0.04, p2 = 12, U(eta) = 0.03*eta, UF(eta) = 0.002*eta];
I used that on the system consisting of ueq2,Teq,eq3,eq4,eq5, where eq4 and eq5 are mentioned in an earlier answer. I had success without problems down to NBT=0.1 when using approxsoln=GUESS together with method=bvp[midrich].
dsolve does all the work including the determination of p2 and phi0.
I have uploaded the worksheet:
MaplePrimes14-07-10odesysInt3.mw

#Since this morning I updated the worksheet to include a loop in NBT from NBT=1 to NBT=0.07. The previous result is used as an approximate solution for the next computation. The start is GUESS.

Try plotting:
plot(log(sec(x)+ tan(x)),x=0..10);
plot(sec(x)+ tan(x),x=0..10,discont=true);

Since sec(x)+tan(x) = (1+sin(x))/cos(x) has the same sign as cos(x) it is negative on the intervals Pi/2..3*Pi/2 and 5*Pi/2..7*Pi/2 etc. Thus the logaritm will have an imaginary part.

These agree:
evalf(Int(log(sec(x)+ tan(x)),x=0..Pi/2));
evalf(int(log(sec(x)+ tan(x)),x=0..Pi/2));
So do these:
evalf(Int(log(sec(x)+ tan(x)),x=3*Pi/2..5*Pi/2));
evalf(int(log(sec(x)+ tan(x)),x=3*Pi/2..5*Pi/2));

But clearly the answer given by Maple on the interval x=0..10 is wrong in the 'int' case.

I suppose you already tried a specific example. Here is one:
pde:=diff(f(t,x),t)=piecewise(x<1,diff(f(t,x),x,x)*f(t,x),diff(f(t,x),x,x));
#debug(`pdsolve/numeric/process_PDEs`);
res:=pdsolve(pde,{f(0,x)=1,f(t,0)=1,f(t,2)=2},numeric);
#You get the error message:
Error, (in pdsolve/numeric/process_PDEs) piecewise can depend only on the independent variables
#which tells the story once you understand that it means that f(t,x) cannot appear inside piecewise at all.
showstat(`pdsolve/numeric/process_PDEs`);
#However, this works:
pde:=diff(f(t,x),t)=piecewise(x<1,sin(x*t),1)*diff(f(t,x),x,x);

You can speed up the procedure by introducing two new dependent variables U and UF defined by
eq4:=diff(U(eta),eta)=u(eta);
eq5:=diff(UF(eta),eta)=u(eta)*phi(eta);
with initial conditions
U(0)=0, UF(0)=0;
Now your system consist of the odes
ueq2,Teq,eq3,eq4,eq5
and the previous boundary conditions together with the two new conditions.
Now define two new local variables Up and UFp by including them in the definition of F0,F1,F2,F3:
Up,UFp,F0,F1,F2,F3:=op(subs(res,[U(eta),UF(eta),u(eta),phi(eta),T(eta),diff(T(eta),eta)]));
The integrals INT0 and INT10 will be Up(1) and UFp(1), respectively.
Otherwise you can procede as before.
I did not, however get any further down than NBT=0.37, when using the continuation
N_bt:=cc*NBT+(1-cc)*1;
together with the fixed values of the other parameters.
The reason I didn't see this first is probably due to the fact (I believe) that originally (a long time ago) your integrals did not have upper limit equal to the endpoint of the interval for the bvp. But now the integral is over 0..1 and so is the bvp. That helps by letting dsolve do all the integration.
Briefly I tried replacing q and LSSolve by including the conditions a[1]=0 and a[2]=0 with boundary conditions at eta=1. dsolve allows undetermined parameters as long as enough boundary conditions are given to determine the parameters (phi0 and p2). That part should be OK. However, I had no immediate success.

The letter 'e' in Maple is nothing but the the letter e.
To get the mathematical constant e use exp(1).
To get e^(-2/7) use exp(-2/7).

Keeping all parameters except NBT as they are in the worksheet and leaving NBT unassigned, I tried with success the following loop:

pt:=[12,0.04]: #Starting initial point
for NBT2 in [1,0.9,0.8,0.7,0.6,0.5,0.45] do
##Using NBT2 as loop variable, leaving NBT unassigned

Q:=proc(pp2,fi0,nbt) option remember; local res,F0,F1,F2,F3,a,INT0,INT10,B;
print(pp2,fi0,nbt);
if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if;
res := dsolve(subs(p2=pp2,phi0=fi0,NBT=nbt,{ueq2,Teq,eq3,u(0)=lambda*D(u)(0),u(1)=-lambda*D(u)(1),D(T)(0)=-1,phi(0)=phi0,(T)(0)=0}), numeric,output=listprocedure,continuation=cc);
F0,F1,F2,F3:=op(subs(res,[u(eta),phi(eta),T(eta),diff(T(eta),eta)]));
INT0:=evalf(Int(F0(eta),eta=0..1));
INT10:=evalf(Int(F0(eta)*F1(eta),eta=0..1));
B:=(-cbf*rhobf+cp*rhop)*INT10+ rhobf*cbf*INT0;
a[1]:=(B-10000*pp2)/10000/pp2; #relative
a[2]:=(INT10/INT0-Phiavg)/Phiavg; #relative
[a[1],a[2]]
end proc:
Q1:=proc(pp2,fi0) Q(pp2,fi0,NBT2)[1] end proc;
Q2:=proc(pp2,fi0) Q(pp2,fi0,NBT2)[2] end proc;
##
RES[NBT2]:=Optimization:-LSSolve([Q1,Q2],initialpoint=pt);
pt:=convert(RES[NBT2][2],list); #To be used at next step
end do;
###An attempt to continue with
pt:=convert(RES[0.45][2],list);
for NBT2 in [0.44,0.43,0.42,0.41,0.40,0.39,0.38] do
etc.
succeeded for NBT2=0.44 only.

#I made a small change after the initial posting: Introducing NBT2 as the loop variable. It was NBT. It didn't seem to matter, though. The danger I saw was in the subs in dsolve, where NBT should be seen as the NBT in the odes.

If the system as given really does have a singularity at about t=0.3300977 then there is nothing you can do about it. It appears from the plot
odeplot(solution,[[t,x(t)],[t,alpha(t)],[t,z(t)],[t,theta(t)]], t=0..0.33009, thickness=2);
that indeed there is a singularity.
There is nothing strange about that for a nonlinear system, just consider as an example the extremely simple equation
ode:=diff(x(t),t)=x(t)^2;
with x(0) = 1:
The solution is found by
dsolve({ode,x(0)=1});
to be x(t) = 1/(1-t), thus has a singularity at t=1. Expressed differently: The maximal solution to this simple ivp has domain of definition (-infinity,1), i.e. not the whole real line.

If for some reason you think your problem should have a solution existing on the whole real line then maybe your mathematical model is wrong.

Incidentally, your 'if theta(t)<>0 then etc.' doesn't make much sense since theta(t) is not given any value by dsolve. Thus theta(t)<>0 is trivially satisfied.

I'm not using the example in your image since I would have to type all that myself.
Not all equations have solutions that can be given in closed form. Thus allvalues will have no choice, but to leave the RootOf alone.
Examples:
solve(x^5+2*x^4+3*x^2+4*x+5=0,x);
solve(x=tan(x),x);
allvalues(%);
#You say your system is linear. Although I really cannot see the system I'm rather confident that your system is not linear. Had it been no RootOf would have appeared.

You can take a look at
?interface
specifically the description of the various names containing 'elision'.

You asked a similar question the other day. My answer to that appears to work for this problem too:
eqs,bcs:=selectremove(has,dsys4,diff);
indets(eqs,function);
eqs[1];
eqs[2];
eqs[3];
eq2:=diff(eqs[2],theta); #Differentiate the equation having diff(h2(theta),theta$4)
isolate(eq2,diff(h3(theta),theta$4));
eq3:=eval(eqs[3],%);
solve({eqs[1],eq2,eq3},{diff(h1(theta),theta$4),diff(h2(theta),theta$6),diff(h3(theta),theta$4)});
params:=indets(eqs,name) minus {theta}=~1;
dsol6 := dsolve(eval({eqs[1],eq2, eq3},params) union bcs,numeric);

You could consider using events:
restart;
with(plots):
odes:=diff(y1(t), t) = y1(t)^2-4*y1(t)+y2(t)*y1(t)-y2(t)+1, diff(y2(t), t) = y1(t), diff(y3(t), t) = -y3(t)+1;
ics:=y1(0) = 0, y2(0) = 0, y3(0) = 0, b(0) = 0;
Sys := {ics, odes};
Sol := dsolve(Sys, numeric, discrete_variables = [b(t)::float], events = [[rhs(odes[2])-rhs(odes[3]), b(t) = t]]);
odeplot(Sol, [t, abs(y2(t)-y3(t))], t = 0 .. 100);
odeplot(Sol, [t, b(t)], 0 .. 100);
odeplot(Sol, [t, abs(y2(t)-y3(t))], t = 65 .. 100);
Sol(5);
Sol(80);

You can get over the problem that convertsys has in isolating the highest derivatives by differentiating the equation having diff(h2(theta),theta$4) as the derivative of highest order:
eqs:=select(has,dsys3,diff);
indets(eqs,function);
eqs[1];
eqs[2]; 
eqs[3];
eq2:=diff(eqs[2],theta); #Differentiate the equation having diff(h2(theta),theta$4) as highest
isolate(eq2,diff(h3(theta),theta$4));
eq3:=eval(eqs[3],%);
#For the new system {eq[1],eq2,eq3} we can solve for the highest derivatives:
solve({eqs[1],eq2,eq3},{diff(h1(theta),theta$4),diff(h2(theta),theta$6),diff(h3(theta),theta$4)});
#We don't have to solve, but did it to see if it can be done!
#Now solving the boundary value problem with all parameters equal to 1 (including `1`):
dsol5 := dsolve(eval({eqs[1],eq2, eq3},{`1`=1,beta=1,kappa=1,chi=1}) union {h1(0) = 0, h1(1) = 0, h2(0) = 0, h2(1) = 0, h3(0) = 1, h3(1) = 1, D(h1)(0) = 0, D(h1)(1) = 0, D(h2)(0) = 0, D(h2)(1) = 0, D(h3)(0) = 0, D(h3)(1) = 0, (D@@2)(h3)(0) = 0,(D@@2)(h3)(1) = 0},numeric);
plots:-odeplot(dsol5,[[theta,200*h1(theta)],[theta,200*h2(theta)],[theta,h3(theta)]],0..1);

An usual name `1` you use, but that is not the problem.
Try solving the odes for the highest derivatives:
eqs:=select(has,dsys3,diff);
solve(eqs,{diff(h1(theta),theta$4),diff(h2(theta),theta$6),diff(h3(theta),theta$4)});
#No solution.
So you have to solve this purely algebraic problem before you can get anywhere.

The problem is not so much about the 'known' option. convertsys cannot isolate theta''(x). In EQT it appears on the left and on the right.
But the known option should be a separate argument to dsolve, i.e. not appear inside the set of equations.
To see that convertsys has problems you can try (with f(x)=1):
dsolve({eval(EQT,f=1),theta(0)=1,D(theta)(0)=.2},numeric,output=listprocedure);

An unrelated problem is that f'''(x) appears on the right in EQT. That should be replaced by the right hand side of EQ before going into the 'known' business. Then f, f', f'' should be replaced by G0,G1,G2, respectively.
#Isolating "by hand":
EQTa:=EQT assuming x<7.5;
EQTb:=isolate(EQTa,diff(theta(x),x,x));
EQTc:=EQT assuming x>7.5;
EQT1:=diff(theta(x),x,x)=piecewise(x<7.5,rhs(EQTb),rhs(EQTc));
EQT1a:=eval(EQT1,EQ);
#remember that G0,G1,G2:=op(subs(subs(res2),[f(x),diff(f(x),x),diff(f(x),x,x)])):
EQT2:=eval(EQT1a,{f(x)=G0(x),diff(f(x),x)=G1(x),diff(f(x),x,x)=G2(x)});
sol:=dsolve({EQT2,theta(0)=1,D(theta)(0)=.2},numeric,known=[G0,G1,G2],output=listprocedure);
plots:-odeplot(sol,[x,theta(x)],0..10);
#In QT you should correct res to res3 and use the known option as above (known=[G0,G1,G2])
fsolve(QT(pp3)=0,pp3=-2..2);
##returns -.3721483873

Using laplace transform on this system fails to solve for x,y,z. You express the laplace transforms of x, y, and z in terms of the laplace transforms of x*y and x(t)/(100+t) etc. Thus you are getting nowhere.
This  nonlinear system cannot be solved by using laplace transforms.
Solve the system numerically.
res:=dsolve({eq1,eq2,eq3,x(0) = 0.5e-1, y(0) = 0, z(0) = 0},numeric);
plots:-odeplot(res,[[t,x(t)],[t,y(t)],[t,z(t)]],0..3000);