I tried using a shooting method for finding f.
I turned the bvp problem for f into an initial value problem with a parameter f2 = f '' (0).
The method then consists in finding the right f2.
As mentioned in my first answer there is a singularity problem because of the coefficient of the highest derivative may become zero. This problem has to be dealt with.
Here I deal with it by changing the aforementioned coefficient when it gets close to zero.
The equation for theta is easily solved after f has been found.
The results seem rather smooth.

MaplePrimes14-01-29b.mw

Version 2: Since it may not be obvious (not to me either) why the approach in the first version is justified, I made another version that takes a somewhat different approach. It has the advantage that the idea is clearer (I hope).
It goes from initially a 3rd order ode to a 2nd order.
The results are, however, not all that different.

MaplePrimes14-01-29b.mw

One problem is that the coefficient of f''' could become zero. Since this is the highest derivative of f occurring in the system you could have a serious problem.
I was not able to find a solution with pi = 0.1, but was able to find a solution for pi as large as 0.0187.
But I admit to being surprised by the fact that the solution f at eta = 10 was clearly greater than the maximum amount allowed for the coefficient of f''' not to become zero!

Note: At the bottom I have shown a different approach.

restart;
Digits:=15:
eq1:=diff(f(eta),eta$3)+1/2*f(eta)*diff(f(eta),eta$2)-pi*(2*f(eta)*diff(f(eta),eta)*diff(f(eta),eta,eta)+f(eta)^2*diff(f(eta),eta,eta,eta)+eta*diff(f(eta),eta)^2*diff(f(eta),eta$2))-K*(diff(f(eta),eta)-1)=0;

eq2:=diff(theta(eta),eta,eta)+1/2*Pr*f(eta)*diff(theta(eta),eta)=0;
bcs:={f(0) = 0, D(f)(0) = 0,D(f)(10) = 1,theta(10) = 1,theta(0)=0};

eq1:=collect(eq1,diff);
#Since f(0) = 0, for this problem not to run into a vanishing coefficient of f''' we need
1-pi*f(eta)^2>0;
#Thus
cd:=abs(f(eta)) < 1/sqrt(pi);
#pi:=0.1;K:=0.2; Pr:=0.7; #No success
pi0:=0.0187; K:=0.2; Pr:=0.7;
eval(cd,pi=pi0);
res:=dsolve({eval(eq1,pi=pi0),eq2} union bcs,numeric);
res(10);
plots:-odeplot(res,[[eta,f(eta)],[eta,theta(eta)]],0..10);

Note. I had success with pi = 0.1 if infinity was replaced by 4.6 instead of 10.
restart;
Digits:=15:
eq1:=diff(f(eta),eta$3)+1/2*f(eta)*diff(f(eta),eta$2)-pi*(2*f(eta)*diff(f(eta),eta)*diff(f(eta),eta,eta)+f(eta)^2*diff(f(eta),eta,eta,eta)+eta*diff(f(eta),eta)^2*diff(f(eta),eta$2))-K*(diff(f(eta),eta)-1)=0;

eq2:=diff(theta(eta),eta,eta)+1/2*Pr*f(eta)*diff(theta(eta),eta)=0;
inf:=4.6:
bcs:={f(0) = 0, D(f)(0) = 0,D(f)(inf) = 1,theta(inf) = 1,theta(0)=0};
eq1:=collect(eq1,diff);
Since f(0) = 0 for this problem not to run into a vanishing coefficient of f''' , we need
1-pi*f(eta)^2>0;
Thus
cd:=abs(f(eta)) < 1/sqrt(pi);
pi0:=0.1; K:=0.2; Pr:=0.7;
eval(cd,pi=pi0);
res:=dsolve({eval(eq1,pi=pi0),eq2} union bcs,numeric,output=listprocedure);
res(inf);
plots:-odeplot(res,[[eta,f(eta)],[eta,theta(eta)]],0..inf);
F,F1,F2,Th,Th1:=op(subs(res,[f(eta),diff(f(eta),eta),diff(f(eta),eta,eta),theta(eta),diff(theta(eta),eta)]));
F(inf);
Fe:=unapply(F(inf)+eta-inf,eta); #Extending f beyond inf
The:=eta->1; #Extending theta beyond inf
Fp:=eta->piecewise(eta<=inf,F(eta),Fe(eta)); # f on 0..infinity
F1p:=eta->piecewise(eta<=inf,F1(eta),1); # f' on 0..infinity
F2p:=eta->piecewise(eta<=inf,F2(eta),0); # f''
Thp:=eta->piecewise(eta<=inf,Th(eta),The(eta)); # theta on 0..infinity
Th1p:=eta->piecewise(eta<=inf,Th1(eta),0); # theta' on 0..infinity
plot([Fp(eta),Thp(eta)],eta=0..6);
plot([F1p(eta),Th1p(eta)],eta=0..6,thickness=3); #theta' not quite continuous at inf.
plot(F2p(eta),eta=0..6,thickness=3);
Clearly the extensions satisfy eq1 and eq2:
eval(lhs~([eq1,eq2]),{f(eta)=Fe(eta),theta(eta)=The});


 

Here is a solution using dsolve with the parameters option.
It seems to me that it solves the problem for N[bt] between 0.2 and 10 and for lambda between 0 and 0.2.

Let me know if you have any questions.

MaplePrimes14-01-23o.mw

restart;
ode := 2*x*y(x)+(x^2+y(x)^2)*(diff(y(x), x));               
dsolve(ode,implicit);
map(t->exp(-3*t),%);
simplify(%);

Assuming now that the model under consideration is Bloch-Grûneisen
http://en.wikipedia.org/wiki/Electrical_resistivity_and_conductivity

with T replaced by T-T0, though, then it can be written

u:=rho[0]+alpha*((T-T0)/theta)^5*(Int(x^5/4/sinh(x/2)^2 ,x = 0 .. theta/(T-T0)));

and can be solved like this (where again I leave out the data copied from this forum):

restart;
Digits:=15:
with(Statistics):
rho[0]:=1.33;
T0:=6;
temp1 and rest1 left out
temp2:=convert(temp1,Vector);
rest2:=convert(rest1,Vector);
plot(temp2,rest2,color=blue); pdata:=%:
T1,T2:=(min,max)(temp2);
(min,max)(rest2);
u:=rho[0]+alpha*((T-T0)/theta)^5*(Int(x^5/4/sinh(x/2)^2 ,x = 0 .. theta/(T-T0)));
f:=proc(T,alpha,theta) evalf(rho[0]+alpha*((T-T0)/theta)^5*Int(x^5/4/sinh(x/2)^2 ,x = 0 .. theta/(T-T0),epsilon=1e-5,method = _NCrule)) end proc;
res2:=NonlinearFit(f,temp2,rest2,output=[parametervalues, residuals],initialvalues=[0.065,4]);
res2[1];
plot(eval(u,{alpha=res2[1][1],theta=res2[1][2]}),T=T1..T2); pmodel:=%:
plots:-display(pmodel,pdata);
## It should be noticed that parameter ranges have been left out. This allows for methods not using the gradient.

MaplePrimes14-01-22N.mw

After some experimentation (in particular by using my favorite tool: animation) I found that an initial point near [15,0.021] ought to work.
In fact I used a shortened version of your system using that phi can be expressed in terms of T. I didn't have to use continuation.
However, the point found by that version [16.7806114174589,0.0241138678090662] ought to work in your original version.
It might help in the procedure Q to make a[2] 500 times larger (scaling some).
If LSSolve ends with "complex value encountered" then you can use the last printout values and check Q.

Incidentally, before you copy and paste code I think you should remove output. In Maple use Edit/Remove Output/From Worksheet.
This makes it easier to try your code.

Assuming now that the missing multiplication sign in
.../(exp(x)-1)(1-exp(-x))
is meant as
.../(exp(x)-1)*(1-exp(-x))
and not as
.../((exp(x)-1)*(1-exp(-x)))
then your integral is easily found by Maple and the result only contains elementary functions.
(The integral is also found with the other interpretation, but involves polylogs.)
Copying and pasting the data as temp1 and rest1 they come out as matrices with one column.
I convert them to vectors.
I have uploaded the full version, but here temp1 rest1 are left out.
To be noticed is that your parameter ranges were not chosen all that well.
restart;
Digits:=15:
with(Statistics):
rho[0]:=1.33;
T0:=6;
As copied:
temp2:=convert(temp1,Vector);
rest2:=convert(rest1,Vector);
plot(temp2,rest2,color=blue); pdata:=%:
T1,T2:=(min,max)(temp2);
(min,max)(rest2);
u:=rho[0]+alpha*(T-T0)*((T-T0)/theta)^5*(int(x^5/(exp(x)-1)*(1-exp(-x)) ,x = 0 .. theta/(T-T0)));
simplify(u):
w:=collect(%,exp,factor);
plots:-animate(plot,[eval(w,alpha=0.07),T=T1..T2],theta=100..500,background=pdata);
res2:=NonlinearFit(w,temp2,rest2,T,output=[parametervalues, residuals],parameterranges=[alpha=0.05..0.08,theta=100..200],initialvalues=[alpha=0.065,theta=450]);
plot(eval(w,res2[1]),T=T1..T2):
plots:-display(%,pdata);

MaplePrimes14-01-22N.mw

The L in NonLinearFit should be l (lower case): NonlinearFit.
Do not use T[0] when T is in use on its own, use T0.
The model function has an integral which Maple cannot do.
You will have to use the operator form: first argument something like this:

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;
#Testing f:
f(7,.0027,50);

The following method has the advantage of applying to equations of the form
diff(x(t),t)=g(t)*x(t-a)+h(t)
where g and h are known functions.
Thus it applies to the equation diff(x(t),t)=x(t-a)+2*t as well as to the equation diff(x(t),t)=16*t*x(t-a) mentioned in one of your remarks.

restart;
ode:=diff(x(t),t)=g(t)*x(t-a)+h(t);
Let H be an antiderivative of h then
eq:=x(t) = x((n-1)*a)+Int(g(tau)*x(tau-a), tau = (n-1)*a .. t)+H(t)-H((n-1)*a);
eq1:=IntegrationTools:-Change(eq,tau=s+a,[s]);
Let X[n](t) = x(t)  for t=(n-1)*a..n*a.
Also taking x((n-1)*a)=X[n]((n-1)*a) = X[n-1]((n-1)*a) for all n we have

eq2:=subs(x(t)=X[n](t),x(s)=X[n-1](s),x((n-1)*a)=X[n-1]((n-1)*a),eq1);
Now suppose x(t) is known on [-a..0]. This means that X[0](t) is known.

X[0]:=t->piecewise(t<0,0,t=0,x0);
h:=t->2*t;  
H:=unapply(int(h(tau),tau=0..t),t);
g:=t->1;
a:=3/10:
N:=4:
for n from 1 to N do
   X[n]:=unapply(simplify(value(rhs(eq2))),t)
end do;
S:=op~([seq([t<i*a,X[i](t)],i=0..N)]):
p:=unapply(piecewise(op(S)),t);
p(t);
eval(p(t),x0=1);
plot(eval(p(t),x0=1),t=-a..1,discont=true,thickness=2);

There is a help page
?Intat

Try
res1:= -Pi+.2+3.0*x+3.2*t+2*Intat(2.25*_f/sqrt(-1.5*_f*(-1.6+4.0500*_f^2+3.7125*_f)), _f = 0);
value(res1);
res2:= -Pi+.2+3.0*x+3.2*t+2*Intat(2.25*_f/sqrt(-1.5*_f*(-1.6+4.0500*_f^2+3.7125*_f)), _f = 1);
value(res2);

The delay differential equation you mentioned must be x'(t)=x(t-3/10)+2*t and not x'(t)=x(t-1)+2*t as you wrote.

Here I find the exact solution to your intended delay ode x'(t)=x(t-3/10)+2*t using laplace transformation.

restart;
with(inttrans);
a:=3/10: #Any a > 0 would work similarly.
ode:=diff(x(t),t)=x(t-a)+2*t;
ode1:=subs(x(t-a)=Heaviside(t-a)*x(t-a),ode); #Assuming x(t) = 0 for t < 0.
dsolve({ode1,x(0)=x0},x(t),method=laplace); #OK, but not finished
#So instead:
laplace(ode1,t,s);
solve(%,{laplace(x(t),t,s)});
L:=subs(x(0)=x0,op(%));
subs(exp(-a*s)-s=-s*(1-exp(-a*s)/s),L);
subs(1/(1-exp(-a*s)/s)=Sum(exp(-n*a*s)/s^n,n=0..infinity),%);
L1:=combine(%);
X:=convert(rhs(invlaplace(L1,s,t)),factorial);
_EnvUseHeavisideAsUnitStep := true;
convert(eval(X,{x0=1,infinity=3,Sum=add}),piecewise,t);
expand(%);
That is your T1!
plot(eval(X,{x0=1,infinity=3,Sum=add}),t=0..1);

Does this work as you intend?
A:=Matrix(3,symbol=a);
B:=Matrix(2,4,symbol=b);
LinearAlgebra:-DiagonalMatrix([A,B]);

evalc works better.
Be aware that in making new assumptions on a name you should use additionally. In your case the assumptions made on mu. However, the assumption mu>0 is irrelevant, so it has no consequence.
Also it helps to make the assumption directly as mu^2 - 4*L*k<0.

restart;
assume(a < 0);
convert(cosh(sqrt(a)), sincos);
evalc(cosh(sqrt(a)));
restart;
assume(K > 0);
#assume(mu > 0);
#additionally(mu^2 - 4*L*k<0);
assume(mu^2 - 4*L*k<0);
assume(t > 0);
convert(cosh((1/2)*t*sqrt(mu^2-4*L*k)/L), sincos); #Doesn't do anything
evalc(%);
about(mu);

This is not possible without using artificial devices. The reason is that even if you imagined having accomplished the fact, i.e. obtained

(x+y)/2;

Then Maple will automatically simplify this to x/2 + y/2 as can be seen by enclosing the expression in unevaluation quotes:

'(x+y)/2';

Questions similar to yours have come up rather often.

By artificial devices I mean something like this:

u:=x/2+y/2;
subs(1/2=``(1/2),u);
factor(%);

Help is available in
?pdsolve,numeric,education
?pdsolve,numeric,methods
Take a look at the examples at the end in ?pdsolve,numeric,education:
You can adjust the very first example to your situation by changing the pde and the ibc to
PDE := diff(u(x,t),t)=diff(u(x,t),x,x);
IBC := {u(x,0)=sin(2*Pi*x), u(0,t)=0,u(2,t)=0};
You need to give the spacestep and timestep explicitly to avoid using the default as in
res1 := pdsolve(PDE, IBC, numeric, time=t, range=0..2,
                 method=Euler,timestep=.1,spacestep=.4); #first case
You have to decide if method Euler is what you want.
Now just continue doing as in the example defining errorfunc (at t=.5) and maxerr (where i=0..200).