IQ:=Vector([seq(Q,Q=-1.5..1.5,.1)]);
DT:=Vector([seq(Dtheta[Q],Q=-1.5..1.5,.1)]);
plot(IQ,DT);

You could perhaps use a spline function for interpolating the data. Here is an example where I construct the data myself.

X:=:
Y:=sin~(X):
plot(X,Y);
S:=CurveFitting:-Spline(X,Y,x);
plot(S,x=0..Pi);
int(S,x=0..Pi);
int(S,x);
plot(%,x=0..Pi);
int(1/(S+1),x=0..Pi);
simplify(%);
int(1/(sin(x)+1),x=0..Pi);
plot(diff(S,x),x=0..Pi);

#If your data was generated in Maple by solving an ode (as in your previous question about writing data to a file) it is probably better to work directly with that ode.

The solutions you found must be used for building the solution to the problem.

Begin by noticing that F0(r,t) =135 satisfies the PDE as well as the boundary condition F0(0.00004,t)=135.
Thus it suffices to find an F1(r,t) satisfying F1(0.00004,t)= 0 as well as the initial condition F1(r,0) = 25-135.
Then your solution is F = F0 + F1.

F1 must consist of infinitely many terms of the form found by Maple. Require each to satisfy the boundary condition. Since BesselY(0,x) -> -infinity as x->0 we exclude contributions from BesselY.
The boundary condition then requires BesselJ(0, 0.00004*lambda) = 0, where I have set c[1]= -lambda^2.
Thus we need the (infinitely many) positive zeros for BesselJ(0,x). The j'th zero can be found in Maple by
BesselJZeros(0, j);
lambda[j] will be that zero divided by r0=0.00004.

Now write F1(r,t) = Sum(c[i]*exp(-lambda[i]^2*k*t)*BesselJ(0,lambda[i]*r),i=1..N), where k=.102*10^(-6) and N is large (10 or 20, in fact it ougth to be infinity).
Determine the c[i]'s by using orthogonality conditions:

int( BesselJ(0,lambda[i]*r)*BesselJ(0,lambda[j]*r)*r, r=0..r0) = 0 for i<>j.

In other words: Maple doesn't give you the solution immediately. You have to work some.

Here is an example using ExportMatrix combined with output = Array. This is probably the easiest way to do it.

sys:=diff(x(t),t)=1-(beta+1)*x(t)+alpha*x(t)^2*y(t), diff(y(t),t)=beta*x(t)-alpha*x(t)^2*y(t);
inits:=x(0)=2,y(0)=0;
alpha:=1: beta:=2.5:
A:=Array([seq(k/10,k=0..200)]):
sol:=dsolve({sys,inits},numeric,output=A);
plots:-odeplot(sol,[x(t),y(t)],0..20);
sol[2,1];
ExportMatrix("F:/MapleDiverse/brussel.dat",sol[2,1]);
restart;
M:=ImportMatrix("F:/MapleDiverse/brussel.dat");
M[1,1..3];
M[-1,1..3];

I don't see any sol2, but I see sol. You are missing a multiplication sign after 0.2. I don't get any error with your code though with or without the multiplication sign.

Edited: Sorry, because I couldn't see the whole text I didn't get your question right. Copying the code and pasting it in a worksheet did give me the whole thing though.

I tried in Maple 12, 15, and 16. No problem. Same answer. I have no access to Maple 14.

If you change to

de1 := dsolve(d1, numeric,output=listprocedure);
de2 := dsolve(d2, numeric,output=listprocedure);
de3 := dsolve(d3, numeric,output=listprocedure);

and

param:={G = 2, B = 2, beta = .2, alpha = 2};
P1 := eval((alpha+3*beta*(diff(f(y), y, y))^2)*(diff(f(y), y, y, y))+G*theta(y)+B*phi(y),param);
P2:=subs(de1,P1);
P3:=evalf(Int(P2,0..1));

then you get a result for P3. But since you set Q at the beginning to 2, it so remains.

And there is no need for the statements with(Student[NumericalAnalysis]); with(PDEtools, casesplit, declare);

This problem is very similar to one I answered on MaplePrimes May 9, 2012.

You can turn the bvp problem into an initial value problem with parameters:

restart;
with(plots):
b := 7;
eq1:=diff(f(eta), eta, eta, eta)-M^2*(diff(f(eta), eta))-(diff(f(eta), eta))^2+m1*f(eta)*(diff(f(eta), eta, eta)) = 0;
eq2:=diff(theta(eta), eta, eta)+m1*Pr*f(eta)*(diff(theta(eta), eta)) = 0;
#Notice that the value of M has been changed from your M=2 to M=1/2 in order to actually get dual solutions, cf. the pdf article you mention.
param:= {m1 = 1, Pr = 0.7, M = 1/2};
res := dsolve({eval(eq1,param), f(0) = S, D(f)(0) = -1, (D@D)(f)(0) = p2}, numeric,output=listprocedure,parameters=[S,p2]);
F0,F1,F2:=op(subs(res,[f(eta),diff(f(eta),eta),diff(f(eta),eta,eta)]));
p:=proc(s,pp2)
res(parameters=[s,pp2]);
F1(b)
end proc;
S0:=1.75:
#curry(p,S0) is the function pp2->p(S0,pp2).
#Plotting to find for which values of f''(0) we have F1(b) = 0, i.e. f'(b)=0:
plot(curry(p,S0),0..3);
plot(curry(p,S0),0.5..1.5);
#We see two solutions:
pm1:=fsolve(curry(p,S0),1);
pm2:=fsolve(curry(p,S0),0.7);
res(parameters=[S0,pm1]);
F1(b);
odeplot(res,[eta,diff(f(eta),eta)],0..b);
PLf1:=%:
resTh:=dsolve({eval(eq2,param union {f(eta)=F0(eta)}), D(theta)(0) = -1, theta(b)=0}, numeric,output=listprocedure,known=F0);
odeplot(resTh,[eta,theta(eta)],0..b);
PLTh1:=%:
res(parameters=[S0,pm2]);
F1(b);
PLf2:=odeplot(res,[eta,diff(f(eta),eta)],0..b,color=blue);
display(PLf1,PLf2);
resTh:=dsolve({eval(eq2,param union {f(eta)=F0(eta)}), D(theta)(0) = -1, theta(b)=0}, numeric,output=listprocedure,known=F0);
odeplot(resTh,[eta,theta(eta)],0..b,color=blue);
PLTh2:=%:
display(PLTh1,PLTh2);
###IMPORTANT: Do the computation of resTh before you change the parameters in res as that will also change F0 (and F1, F2). 

Your functions q and p are well defined on the real line:

lambda/mu*q(-1);
lambda/mu*q(-2);
lambda/mu*q(1.2345);

If you want to restrict its definition to the nonnegative integers you can do:

q := (n::nonnegint) -> (1-r[0])*r[0]^n;

Since there is success for da4 in finding a solution, you can use that to find solutions for the other values of S:

da3:=dsolve(d3, numeric,approxsoln=da4);
da2:=dsolve(d2, numeric,approxsoln=da3);
da1:=dsolve(d1, numeric,approxsoln=da2);

This gives you one solution for f and theta for your values of S. So that is a start.

But what makes you think that there are two solutions for each set of values for the parameters?

The answer may not tell you much, but it is not absurd.

Quoting from the help page for dsolve,details:

"For linear ODEs, when dsolve is not able to find a solution or a reduction of order, an answer using DESol is returned. DESol structures can also appear in answers involving reductions of orders of linear ODEs. Although DESol structures do not contain more information than the ODE itself, they can be useful in further computations since Maple can manipulate them, for example, by expanding in series, or by simplifying. See DESol."

s:=CurveFitting:-Spline(Cl_alpha,x):
S:=unapply(s,x):
ynew:=S~(xnew);
ynew[5,1];
ynew(5);

Here is an example. I make up my own data:
X:=[i$i=0..30];
Y:=evalf[2](sin~(X));
s:=CurveFitting:-Spline(X,Y,x):
plot([s,sin(x)],x=0..30);
S:=unapply(s,x):
S(12.3456);
#Alternatively:
eval(s,x=12.3456);

#I noticed that this is closely related to your previous question, which I saw after answering this.

You could open the .mla file in a worksheet. Maple then uses 'march'. After that use eval(myproc);

Or you could in a worksheet do something like this:

march('open',"F:\\Mat\\test.mla");
march(list,"F:\\Mat\\test.mla");
showstat(myproc);
eval(myproc);

This does indeed appear to be difficult.
I didn't have any success with the direct approach, so tried dsolve on an initial value problem with y2=y''(0) and y3=y'''(0) as parameters. There I had some success when Digits = 10 and Digits = 12. However, I didn't have success with Digits=15. So is what you will see here a solution, or is it not?

restart;
Digits:=12:
deq := diff(y(x), x$4) = -.364*diff(y(x), x)^2*diff(y(x),x,x)+k;
bc := y(0) = 0, D(y)(0) = 0, y(L) = 0, D(y)(L) = 0;
k0:=1.74*10^(-7);
L:=1000;

ics:=y(0) = 0, D(y)(0) = 0, (D@D)(y)(0)=y2,(D@@3)(y)(0)=y3;
dsolpar := dsolve({ics, eval(deq,k=k0)}, numeric,parameters=[y2,y3],output=listprocedure,maxfun=0,range=0..L);

Y,Y1:=op(subs(dsolpar,[y(x),diff(y(x),x)]));
p1:=proc(y2,y3)
dsolpar(parameters=[y2,y3]);
Y(L)
end proc;
p2:=proc(y2,y3)
dsolpar(parameters=[y2,y3]);
Y1(L)
end proc;
fsolve([p1,p2],[0.005..0.015,-0.0002..-0.0001]);
#[0.7898528584778143356e-2, -0.1372579321779395876e-3]
dsolpar(parameters);
Y(L);
Y1(L);
with(plots):
display(Array([odeplot(dsolpar,[x,y(x)],0..L,refine=1),odeplot(dsolpar,[x,diff(y(x),x)],0..L,refine=1)]));
display(Array([odeplot(dsolpar,[x,y(x)],0.9*L..L,refine=1),odeplot(dsolpar,[x,diff(y(x),x)],0.9*L..L,refine=1)]));

I actually began with the direct approach which worked fine with L = 250 and k = k0 or equivalently with the constant k= k0/256 and L=1000. Also L=1000 and k up to k0/176 worked.