@tobbie247 You can try dsolve as in the following, where I had success with infinity = inf = 28, but not with inf = 29.
For my convenience I used x instead of eta.

ode1:=diff(f(x),x$3)+3*f(x)*diff(f(x),x,x)-2*diff(f(x),x)^2+theta(x)-m*diff(f(x),x)=0;
ode2:=diff(theta(x),x,x)+3*Pr*f(x)*diff(theta(x),x)+s*theta(x)=0;
params:={Pr=1,s=1,m=2};
bcs:={f(0)=0,D(f)(0)=0,theta(0)=1,D(f)(inf)=1,theta(inf)=0};
sys:=eval({ode1,ode2},params);
inf:=28:
res:=dsolve(sys union bcs,numeric);
plots:-odeplot(res,[[x,diff(f(x),x)],[x,theta(x)]],0..inf);
plots:-odeplot(res,[[x,f(x)],[x,theta(x)]],0..inf);

##What is the meaning of A and B in the discussion above in relation to the system of odes?

@rlewis Did you not just copy and paste all the lines I gave?

I just did that myself directly from my answer above. It ran without problems in Maple 2015.1.

Which Maple version are you using?

## Note: It seems that you are using a version behaving like Maple 12, see below. ###############

##########
I just tried in Maple 15. There the third ball is white (!!), thus cannot be seen!
You can change that by setting the color, either a common color as in color=blue or 3 different colors as shown in the revised lines below. There I also add the option axes=box. In Maple 2015 this is the default, in Maple 15 the default is axes = none.

restart;
with(plots):
f1:=cos(t)/4: f2:=sin(t)/5: f3:=sin(2*t+1)/6:
pts:=[[0,0,f1],[1+f2,1,1],[-1,0,f3]];
ptsL:=[pts[],pts[1]];
pointplot3d(eval(pts,t=0),symbol=solidsphere,symbolsize=50,color=[blue,red,green]); p1:=%:
pointplot3d(eval(ptsL,t=0),connect,thickness=3,color=black); p2:=%:
display(p1,p2);
#Now animating from t = 0 to t = 50:
animate(display,['pointplot3d'(pts,symbol=solidsphere,symbolsize=50,color=[blue,red,green],axes=box),'pointplot3d'(ptsL,connect,thickness=3,color=black)],t=0..50,frames=100);

###################### Maple 12 #######
In Maple 12 I see the error message you report. Instead of the previous method of animation you can do:
restart;
with(plots):
f1:=cos(t)/4: f2:=sin(t)/5: f3:=sin(2*t+1)/6:
pts:=[[0,0,f1],[1+f2,1,1],[-1,0,f3]];
ptsL:=[pts[],pts[1]];
pointplot3d(eval(pts,t=0),symbol=solidsphere,symbolsize=50,color=[blue,red,green]); p1:=%:
pointplot3d(eval(ptsL,t=0),connect,thickness=3,color=black); p2:=%:
display(p1,p2);
#Now animating from t = 0 to t = 50:
p:=tt->display(pointplot3d(eval(pts,t=tt),symbol=solidsphere,symbolsize=50,color=[blue,red,green],axes=box),pointplot3d(eval(ptsL,t=tt),connect,thickness=3,color=black));
animate(p,[t],t=0..50,frames=100);

@Kanellopoulos There are a lot of things to vary in these kinds of numerical problems, so I wouldn't dare conclude much else than for the settings we are dealing with in Carl's worksheet, results don't seem reliable.
I tried (using sysS as described above) to look at relative residuals. By that I mean the difference between the left and right hand sides of the pdes divided by the sum of their absolute values.
Since sysS is a set as I defined it, I replaced it by
sysS1:=convert(sysS,list);
so that the order is known.
Using the relative residuals meant replacing (lhs-rhs) by
((lhs-rhs)/(abs@rhs+abs@lhs))
in the definition of Residuals.
The results were rather discouraging. I found that in 8 pairs of residuals one of them (always the second: pde2) the relative residual had absolute value 1.

@Carl Love I tried your worksheet with the following change.
Using your sys I did
sysS:=solve(value(sys),{diff(h(x,t),t),diff(u(x,t),t)});

Then I used that instead of sys in the rest of the worksheet, i.e. in defining sol and in defining Residuals.

I got residuals over 0.1 in absolute value in 7 of the 18 (X,T) cases. One even over 3.

@tomleslie I tried making a new document in document mode and 2d input, just 3 lines. I copied this and pasted it into a fresh worksheet (worksheet mode and 1d input).
There I converted the pasted code into 1d-input. That worked, but there are no prompts or execution group delimiters and the area covered by the lines seems to be a document block.
This can be changed by using the menu item Format/Remove Document Block.
Unfortunately this meant that (in my case) 3 lines of
print(??); # input placeholder
also appeared.
#############
I just had a look at the discussion about converting a file made in document mode and 2d input into a file in worksheet mode and 1d input.
http://www.mapleprimes.com/questions/203080-How-To-Convert-Document-In-Document

I had occasion to do this the other day
http://www.mapleprimes.com/questions/204487-Too-Lengthy-Parametric-Equations

It would have been easier to do the conversion by exporting to .mpl as suggested by Alejandro Jakubi in the first reference.

@cuongtd The output from allvalues(solve....) is not a numerical result, but is the exact result.
As it turns out all 32 solutions (16 real and 16 imaginary) can be expressed in terms of the 16 roots of a polynomial of degree 16. Once you use evalf on that exact result numerical rootfinding is used, but since the problem just lies with the polynomial surely all roots are found.
To see the result before using evalf just change the colon in the first line below into a semicolon: It is huge!
#I'm assuming that eqs and eqs2 are defined as I did in my answer. Still using Digits=30.

resA:=allvalues([solve(eqs union eqs2)]): ##Sequence of length 16 each consisting of a list of two solutions
nops(resA[1]); # A pair of 2 solutions
nops([resA]); #16 pairs
nops(ListTools:-Flatten([resA])); #32 solutions
indets(resA[1],specfunc(RootOf)); #Notice the two arguments to RootOf: a polynomial of degree 16 and an index
ROF:=indets([resA],specfunc(RootOf)); #All of the RootOf's
nops(ROF); #16
pol:=map2(op,1,ROF); #The same polynomial for all
map2(op,2,ROF); #The 16 indices
#Now using evalf as I originally did.
res:=evalf(resA):
resR:=remove(has,[res],I); #Looking at the real results only
##Comment: It so happens that each two in the 16 pairs are either both real or both imaginary. I checked that.
map(nops,resR); #8 pairs
resRflat:=ListTools:-Flatten(resR);
nops(%); #16
##Graphical illustration: Animation showing one solution at a time.
##The 4 points in each picture are (ci,si), i=1..4.
plts:=seq(plots:-pointplot(subs(resRflat[j],[seq([cat(c,i),cat(s,i)],i=1..4)]),color=[red,blue,green,maroon],symbolsize=25),j=1..nops(resRflat)):
plots:-display(plts,insequence);
##Note: Recall that my (c4,s4) is the point (cos(t4),sin(t4)), where t4 = t+t1.

@mapleus I notice that A(0) is imaginary. Is that intended?

eval(A(0),r(0)=0.3-0.01);
#Outout: 9.330399019*I

Before trying a loop like that it makes sense to make it work for just one step!
That is, make dsolve come out with some result. I got an error.

Besides that I see one problem right away:
subs(F,r(s)) is itself a procedure, so you should just have
r_ravn:=subs(F,r(s));

@nafiqah38 As you can see yourself there is nothing at all in your reply.

Give us your system with boundary conditions etc. either as pasted code (text) or as an uploaded worksheet. For the latter use the thick green arrow in the editor.

@FMEHR The approach I showed above works with an appropriate extra condition.

restart;
dsys3 := {-901.7765286*(diff(f1(x), x, x))+3.153274902*10^(-10)*(diff(f2(x), x, x, x))+0.4807756008e-5*(diff(f2(x), x))-27.92641797*(diff(f3(x), x))+18074.25056*f1(x)-3.100642108*(diff(f3(x), x, x))*(diff(f3(x), x))+7.74301535*(diff(f3(x), x))*f3(x) = 0, 0.1873004296e-2*(diff(f3(x), x, x, x, x))+(diff(f3(x), x, x))*(-0.1142850079e-2-0.1587291776e-1*P)-.1797931318*(diff(f3(x), x, x))+13.96314014*(diff(f1(x), x))+1.822500000*10^(-10)*(diff(f2(x), x, x))-4.888800000*10^(-7)*f2(x)+22.72337812*f3(x) = 0, 0.2249991345e-1*(diff(f2(x), x, x, x, x))-227.1861332*(diff(f2(x), x, x))+29537.78948*f2(x)-2.841067252*10^(-11)*(diff(f1(x), x, x, x))-3.39267254*10^(-7)*(diff(f1(x), x))+7.650000000*10^(-11)*(diff(f3(x), x, x))-1.400000000*10^(-7)*f3(x)-5.500000000*10^(-7)*f3(x)^2+9.278333333*10^(-13)*(diff(f3(x), x))^2+3.300000000*10^(-13)*(diff(f3(x), x, x))*f3(x) = 0, f1(0) = 0, f1(1) = 0, f2(0) = 0, f2(1) = 0, f3(0) = 0, f3(1) = 0, (D(f2))(0) = 0, (D(f2))(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0};
 
indets(dsys3,name);
for i to 3 do dsys3[i] end do;
diff(dsys3[2],x);
eq:=subs(solve(dsys3[3],{diff(f2(x),x$4)}),%);
sys:=solve({dsys3[1],dsys3[3],eq},{diff(f1(x),x$3),diff(f2(x),x$4),diff(f3(x),x$4)});
dsys3[4..-1];
cdn:=eval(convert(dsys3[2],D),x=0);
#Using the extra condition f3'''(0)=-1 we have success with
res:=dsolve(sys union dsys3[4..-1] union {cdn, (D@@3)(f3)(0)=-1},numeric);
res(0); # Eigenvalue P = -42.85...
#If we plot the three together two of them are drowned by the largest (f3). So scale:
plots:-odeplot(res,[[x,100*f1(x)],[x,10^11*f2(x)],[x,f3(x)]],0..1,thickness=3,legend=[400*f1, "10^11*f2", f3],size=[1024,600]);
#Or in 3 separate graphs
use plots in
  display(Array([seq(odeplot(res,[x,cat(f,i)(x)],0..1),i=1..3)]))
end use;
## Finding other P's
## You can try changing 'a' in the extra condition (D@@3)(f3)(0)= a (a<>0). I tried a = 0.4 as an example (and started with Digits:=30). In the linear case this change should not be expected to change matters a priori since if f = (f1,f2,f3) is an eigenvector (function) so is a scalar multiple of f. But your system is nonlinear.
I got P = -41.10... and the graphs were different.

Take a look at the thread:

http://mapleprimes.com/questions/204463-Nonlinear-Ode-Equations

where someone else asked a very similar question just using different constants, but the very same notation: dsys3, f1,f2,f3.
My guess is that it is a fellow student. Talk to him.

@emmantop What I meant is illustrated in this very simple example:
restart;
ode:=diff(x(t),t,t,t)+x(t)/(3-t)=0;
dsolve({ode,x(0)=0,D(x)(0)=1,x(3)=0},numeric); #Error because of the denominator 3-t.
res:=dsolve({ode,x(0)=0,D(x)(0)=1,x(3)=0},numeric,method=bvp[midrich]);
plots:-odeplot(res,[t,x(t)],0..3);
res(3); #Here we get the values of x, x', x'' at 3
res(0); #Here we get the values of x, x', x'' at 0

@FMEHR The image above was produced by my own unsophisticated solver. I have not yet been able to get that by dsolve.

After observing that it produced a result with f2 = 0 I tried using that directly as follows (assuming that sys and cdn are defined as above):

sys20:=eval(sys,f2(x)=0);
for i to 3 do sys20[i] end do; #Viewing the result
cdn20:=eval(cdn,f2=0);
bcs20:=remove(has,dsys3[4..-1],f2);
#Leaving out sys20[1]:
res0:=dsolve(sys20[2..3] union bcs20 union {cdn20},numeric,approxsoln=[f1(x)=-sin(Pi*x),f3(x)=100*sin(2*Pi*x)]);
res0(0.5);
plots:-odeplot(res0,[[x,f1(x)],[x,f3(x)]],0..1); #Pretty small: Reliable???
#Checking if sys20[1] happens to be satisfied:
plots:-odeplot(res0,[x,lhs(sys20[1])],0..1,thickness=5); #It seems so

@FMEHR Do you have any reason to believe that a nontrivial solution exists?

Possibly something like the following, where f2 is the blue graph, i.e. f2(x) = 0 (almost at least)?

I don't trust the solution above since I cannot make dsolve come out with that answer even if I feed it as an approximate solution. It still returns the trivial solution.