eq1:=diff(f(eta),eta,eta,eta)+f(eta)*diff(f(eta),eta,eta)+beta*(1-diff(f(eta),eta)^2)=0;
eq2:=diff(theta(eta), eta, eta)+pr*f(eta)*(diff(theta(eta), eta))+pr*Ec*(diff(f(eta), eta, eta))^2
+(2-beta)*pr*Q*theta(eta) = 0;
bc:=theta(9)=0,D(theta)(0)=-Nc*(1-theta(0)),f(0)=1,D(f)(0)=lambda,D(f)(9)=1;
Para:={beta=1,Ec=0.5,Q=-0.5,Nc=5,pr=3};

#The approach is straightforward in the sense that it just hands dsolve the boundary value problem: We don't do the shooting or whatever method is used. The 'parameters' option is only available for IVPs, so we make a procedure, which for given lambda returns -D(theta)(0):

mTp0:=proc(lambda1) local res,T1;
   if not type(lambda1,numeric) then return 'procname(_passed)' end if;
   res:=dsolve(eval({eq1,eq2,bc},Para union {lambda=lambda1}),numeric,output=listprocedure);
   T1:=subs(res,diff(theta(eta),eta));
   -T1(0)
end proc;
plot(mTp0(lambda),lambda=-2..2);

But experimenting with different ranges for lambda you will notice jumps from one to another of the multiple solutions we discussed in your earlier postings.

Because the numbers are so huge, take the logarithm elementwise of the vector:

L:= Fit(La-b/t, T, ln~(Kpco2F), t);

A:=<1,3,5>;
B:=10^~A;

There is an obvious solution f(eta) = 1 for all eta and for n > 0.

If you insist on attempting dsolve/numeric you could replace D(f)(12) = 0 with D(f)(0)=0.0000001 or some such number.

Added:

After solving for f ''' (which must be done by dsolve) we can write the equation in the form

diff(f(eta), eta, eta, eta) = -f(eta)*signum(diff(f(eta), eta, eta))*abs(diff(f(eta), eta, eta))^(2-n)/(n+1);

Here we are able to handle n <=1.

Here are two similar ways. I admit that the first one makes me uneasy because of the apparent recursive definition. There is no such worry about the second.
restart;
a:=()->a;
f := a+D(x);
f(t);
convert(f(t), diff);

restart;
a1:=()->a;
f := a1+D(x);
f(t);

I am guessing that you are trying to determine for which values of n and t there are nontrivial solutions.

If so try (edited once again!!):

restart;
ode:=t*diff(phi(x),x$2)-x*diff(phi(x),x)+n*phi(x)=0;
sol := dsolve(ode);
odetest(sol,ode);
eq1:=limit(rhs(sol),x=0,right)=0;
eq2:=eval(rhs(sol),x=1)=0;

LinearAlgebra:-GenerateMatrix([eq1,eq2],[_C1,_C2]);
CharEqn:=LinearAlgebra:-Determinant(%[1])=0;

For concrete values of m you can use 'residue'. Otherwise use the formula

Res(f,z0) = limit(  1/(m-1)!*diff((z-z0)^m*f,z$(m-1)),   z=z0  );
###########
f:=exp(a*z)/z/(z+b)^m;
residue(eval(f,m=2),z=-b);
series(eval(f,m=2),z=-b,3);
1/(m-1)!*diff((z+b)^m*f,z$(m-1));
res:=limit(%,z=-b);
eval(res,m=1);
eval(res,m=2);

I have uploaded a worksheet producing the plot. It is a revised version of the one I included in a comment earlier. The plot was produced using implicitplot several times. The bending over showing 4 solutions is not seen in the paper referred to in your comment. I didn't use the shooting library, but the method is of that type.

MaplePrimes12-09-06B.mw

You could use Analytic in the RootFinding Package:

eq:=x-a*(1-x)* exp(b*x)+c=0;
eq1:=eval(eq,{a=11.85646836, b=1/2, c=9.85646836});
RootFinding:-Analytic(eq1,x,-10-10*I..10+10*I);
eq2:=eval(eq,{a=11, b=1/2*I, c=10});
RootFinding:-Analytic(eq1,x,-10-10*I..10+10*I);

If the dependent variable in your 4th order ode is called x then presumably the second column contains values of x(t) for the values of t in the first column and the third column contains values of x'(t).

As suggested by me and strongly emphasized by acer in another question you posed:
http://www.mapleprimes.com/questions/136957-IntegratingDifferentiating-Numerical-Data

it is better to work with the 4th order ode together with your new system.

Even if you prefer to solve the 4th order ode independently it will be better and easier to keep working with the ode directly using output=listprocedure and using the output procedures as input to your system.

Here is an ad hoc example. I first present the latter approach, i.e. solving the 4th order ode before attacking the system.

eq1:=diff(x(t),t$4)=sin(x(t))-cos(diff(x(t),t,t));
icx:=x(0)=1,D(x)(0)=1,(D@@2)(x)(0)=1,(D@@3)(x)(0)=0;
resx:=dsolve({eq1,icx},numeric,output=listprocedure);
X,X1,X2,X3:=op(subs(resx,[seq(diff(x(t),[t$k]),k=0..3)]));
plots:-odeplot(resx,[seq([t,diff(x(t),[t$k])],k=0..3)],0..10);
#And here is your system (just one equation, but that is irrelevant):
eq2:=diff(y(t),t,t)=-y(t)*X(t)+tanh(X1(t)/X(t));
icy:=y(0)=0,D(y)(0)=1;
#Use the 'known' option:
resy:=dsolve({eq2,icy},numeric,known=[X,X1]);
plots:-odeplot(resy,[t,y(t)],0..10);
#Now the easier approach:
eq2a:=subs({X(t)=x(t),X1(t)=diff(x(t),t)},eq2);
resxy:=dsolve({eq1,eq2a,icx,icy},numeric);
plots:-odeplot(resxy,[t,y(t)],0..10);

#If for some reason you really insist on using the dat-file, you could use spline approximation to get to something like X and X1 above.

You could take a look at pade in the numapprox package. 

Example:

with(numapprox);
pade(tan(x), x=0, [5,6]);

Try this:

restart;

with(plots):
eqs:=diff(f(eta),eta,eta,eta)+f(eta)*diff(f(eta),eta,eta)+(1-diff(f(eta),eta)^2)+A*((1-diff(f(eta),eta))-1/2*eta*diff(f(eta),eta,eta))=0,1/pr*diff(theta(eta),eta,eta)+f(eta)*diff(theta(eta),eta)-theta(eta)*(diff(f(eta),eta))-A*(theta(eta)+eta/2*diff(f(eta),eta))=0;

respar:=dsolve({eqs[1],f(0)=f0,D(f)(0)=f1,(D@D)(f)(0)=f2},numeric,output=listprocedure,parameters=[A,f0,f1,f2]);
F,F1,F2:=op(subs(respar,[f(eta),diff(f(eta),eta),diff(f(eta),eta,eta)]));
p:=proc(N,A,f0,f1,f2)
respar(parameters=[A,f0,f1,f2]);
F1(N)-1
end proc;
Taking infinity=8:
plot(tanh@curry(p,8,.01,-.1,-1.18),0.7..0.9);
pm1:=fsolve(curry(p,8,.01,-.1,-1.18),0.74..0.75);
pm2:=fsolve(curry(p,8,.01,-.1,-1.18),0.76..0.79);
respar(parameters=[.01,-.1,-1.18,pm1]):
pf1:=plots:-odeplot(respar,[eta,diff(f(eta),eta)],0..8):
dsolve({eval(eqs[2],{A=.01,pr=0.7,f(eta)=F(eta),diff(f(eta),eta)=F1(eta)}),theta(0)=1,theta(8)=0},theta(eta),numeric,known=[F,F1]):
pt1:=odeplot(%,[eta,theta(eta)],0..8):

display(pf1,pt1,labels=[eta,typeset(diff(f(eta),eta)," and ",theta(eta))]);
respar(parameters=[.01,-.1,-1.18,pm2]):
pf2:=odeplot(respar,[eta,diff(f(eta),eta)],0..8,color=blue):
dsolve({eval(eqs[2],{A=.01,pr=0.7,f(eta)=F(eta),diff(f(eta),eta)=F1(eta)}),theta(0)=1,theta(8)=0},theta(eta),numeric,known=[F,F1]):
pt2:=odeplot(%,[eta,theta(eta)],0..8,color=blue):

plots:-display(pf2,pt2,labels=[eta,typeset(diff(f(eta),eta)," and ",theta(eta))]);
display(Array([display(pf1,pf2),display(pt1,pt2)]));

By substituting the solution into your equations it is seen that the solution given by Maple is correct.

I suppose you have reason to believe that there are other solutions as well?

(It would have been more helpful for us if you had suppressed the output in the assignments!)

The second element in the output doesn't appear to contain an equation but an algebraic expression expr, but as is often the case in Maple it is the implied equation expr=0. A much simpler example:

sol:=eliminate( {x^2+y^2-1=0, x+y=0}, x);
results in

Compare with

solve({x^2+y^2-1=0, x+y=0}, {x,y});

Maybe you could provide the data, since those seem to be the problem. After all the following, where I construct my own data, works:

restart;
f:=x->A*(1-exp(B*x^C));
A:=2: B:=-0.5: C:=.3:
X:=<seq(0.1*k,k=0..10)>;
Y1:=f~(X);
Statistics[Fit](a*(1-exp(b*x^c)),X,Y1,x,initialvalues=[a=1.5,b=1,c=.5]);
Statistics[Fit](a*(1-exp(b*x^c)),X,Y1,x,initialvalues=[a=1,b=1,c=.8]);
Y2:=Y1+Vector(11, RandomTools:-Generate(float(range=-.1..0.1), makeproc=true));
Statistics[Fit](a*(1-exp(b*x^c)),X,Y2,x,initialvalues=[a=1.7,b=1,c=.5],iterationlimit=100000);
plot(%,x=0..1);
plots:-display(%,plot(X,Y2));