You can use odeplot on an ode system that has the equation as an (implicitly given) solution.
I shall use the example provided by Kitonum.

restart;
ode1:=diff(x(t),t)=D[2](f)(x(t),y(t));
ode2:=diff(y(t),t)=-D[1](f)(x(t),y(t));
diff(f(x(t),y(t)),t);
eval(%,{ode1,ode2}); # 0, so f(x(t),y(t)) = C is a solution.
eq:=x^5+y^5-x*y=1; #Example equation
f:=unapply((lhs-rhs)(eq),x,y);
eval(eq,x=0); #Used to find suitable initial conditions.
res:=dsolve({ode1,ode2,x(0)=0,y(0)=1},numeric);
plots:-odeplot(res,[x(t),y(t)],-0.3..0.64);
plots:-odeplot(res,[x(t),y(t)],-0.3..0.64,frames=50);

I only looked at the asymptotic expansion given in equation (45) in the pdf document.
So all I do here is to use odeplot for graphical evidence of the expansion:

restart;
Digits := 15;
inf := 4;
equ1 := diff(f[0](eta), eta$3)+theta[0](eta);
equ2 := diff(theta[0](eta), eta,eta)+3*f[0](eta)*diff(theta[0](eta), eta);
Bcs1 := f[0](0) = 0, D(f[0])(0) = 0, theta[0](0) = 1, theta[0](inf) = 0, (D@@2)(f[0])(inf) = 0;
##
S1 := dsolve({Bcs1, equ1, equ2}, numeric);
##
plots:-odeplot(S1,[[eta,f[0](eta)-0.5106804*eta+0.261009]],
            caption=typeset(f(eta)-(0.5106804*eta-0.261009)));

Your last line has 10^(-21) , where the corresponding element in Sol has 10^(-24).

Try just doing
select(type,[Sol]*10^19,positive);
sqrt~(%);
## I get [.983503318731514, 0.940901039630522e-2] with Digits at 15, and not nuch different at Digits=10, [0.940901039700659e-2, .983503318863600].
#### General comment:
Solutions for omega2 are heavily dependent on the value of b because of the term
0.8325000000e-4*omega2, which has no g1 in it, just as in this example:

ode:=diff(y(x),x,x)+omega2*1.2345=0;
b:=1:
res:=dsolve({ode,y(0)=0,y(1)=0,D(y)(1)=b},numeric);
res(0);
b:=10:
res:=dsolve({ode,y(0)=0,y(1)=0,D(y)(1)=b},numeric);
res(0);

It is never a good idea to use a name (like theta) and an indexed version of the same (like theta[p] ) in the same session as independent names, as you are doing.
Tom Leslie didn't notice it I'm sure because it is handled fine in Maple 2016, but not in the version you are using, Maple 12. (The problem exists in all versions up to and including Maple 2015).
The remedy is simple: Replace theta[p] by another name like thetap or whatever, just not an indexed version of the other ones.
The title of your question led me to guess the reason. I have checked in Maple 12.
You still have to follow the directions given by Tom Leslie too.
##
## Note 1: The problem is with `dsolve/numeric/bvp/convertsys` in Maple 12.
## There is no such problem with DEtools[convertsys] in Maple 12. That procedure is used for the same purpose of converting to an explicit first order system but for initial value problems.
## Note 2: There should be no problem of using only indexed versions of a name, as e.g. theta[0] and theta[1]. It is the mixture of indexed and unindexed that is a problem.

In Maple 8 you can go to the help page

?plot3d, colorfunc

and then experiment with a color option like this one
color=[cos(5*u),cos(5*u)*cos(2*v),cos(15*u)+0.02*cos(25*v)]
or any other list of 3 functions of u and v.
Alternatively, you can use the simpler color option:
color=cos(5*u)
where the right hand side can be a function of u and v.

You need to give mu a value.

Depending on that value you may have to use the option maxfun=0 (which really means infinity) or maxfun = big value in dsolve/numeric.
Furthermore you may want to give the option numpoints=10000 or (some such high number) to odeplot and may also have to restrict the range -300..300 to something smaller.

You have 3 initial conditions for two second order odes. You need 4 initial conditions.
Try this (no numerics):
dsolve({ode1,ode2,ics});
You will notice that the answers represent infinitely many solutions.
As a numerical example of one of these try

res:=dsolve({ode1,ode2,f(0)=0.3,g(0)=0.7,D(f)(0) =0,D(g)(0) =1},numeric);
plots:-odeplot(res,[[x,f(x)],[x,g(x)]],0..4,view=0..1.5);

I'm pretty sure that the 3D plot you link to is not built on the 2D plots. Those 2D plots are just there to illustrate the vertical shape, the horizontal shape, and the vertical perturbation.
Surely the plot is built as the following.
 

V:=A*exp(-b*abs(r)^1.5)*abs(r)^a;
H:=R0+R1*abs(sin(c*theta/2))^d;
Vp:=1+r^2*A1*sin(p*theta);
params:=[a=0.5,b=0.15,A=2,c=5,d=0.5,R0=0.5,R1=0.5,p=15,A1=0.05];
plot3d(eval([H*r*cos(theta),H*r*sin(theta),V+Vp],params),theta=0..2*Pi,r=-1..1,style=patchnogrid,colorscheme=["zgradient", ["Red", "Blue"]],lightmodel="light4");
##
## I now notice that you are using Maple 8. Then you have to remove two options in plot3d. 
##Use just

plot3d(eval([H*r*cos(theta),H*r*sin(theta),V+Vp],params),theta=0..2*Pi,r=-1..1,style=patchnogrid);

##and add whatever you like in Maple 8

There is color work to be done. Furthermore surely this could be done using Explore.
The plot is this:

with(StringTools):
WildcardMatch( "?e?d?", "sends" );
L:=["abcde","fghij","fends","sends"];
select(s->WildcardMatch( "?e?d?", s ),L); #Not Select from StringTools

When supplying necessary information: range and time variable:
PDES := pdsolve({PDE1}, {BC1, IBC1}, numeric, timestep = .1, spacestep = .1, time = t, range = 0 .. 1);
I get the error
Error, (in pdsolve/numeric/par_hyp) input system is too far from a 'standard' form (see ?pdsolve,numeric for more detail)

A little debugging in `pdsolve/numeric`, `pdsolve/numeric/process_PDEs`, and `pdsolve/numeric/par_hyp` shows that the error comes from the latter procedure's lines 88-91. The INFO table assignments used here are done in `pdsolve/numeric/process_PDEs`.
I tried changing  the equation
eta*(diff(R(x, t), t))+diff(U(x, t), x)-R(x, t) = 0
into
eta*(diff(R(x, t), t))+diff(U(x, t), x)-diff(R(x, t),x) = 0

and left the rest as they were. Then everything was just fine.

You have an ode, not a pde.

ode:=(diff(y(x), x)+1)*exp(y(x)) = 4;
dsolve({ode,y(1)=0});

You could try using the result from dsolve/numeric/bvp to find the initial values at one of the two endpoints, in your case they are -0.5 and 0. Indeed you could use any point inside the interval instead. Then use these as initial conditions in dsolve/numeric.
I have tried this in an extremely simple case:

ode:=diff(y(x),x,x)+y(x)=0;
res:=dsolve({ode,y(0)=0,y(Pi/2)=1},numeric); #This is the bvp solution
plots:-odeplot(res,[x,y(x)-sin(x)],0..Pi/2); #Just checking how good it is.
ics:=eval(convert(res(0),D),x=0)[2..-1]; #Will be taken as initial conditions
res_ics:=dsolve({ode,op(ics)},numeric); #NOT using dsolve/numeric/bvp
plots:-odeplot(res_ics,[x,y(x)-sin(x)],-2*Pi..3*Pi); #Not too bad

I find it a bad idea to use lower case l, since it looks quite like 1.

You have

(D(f))(0) = l+b*(D@@2)(f)(0);

where l := [1, 2, 3], which clearly is a problem.

I don't know why, but dsolve/numeric/bvp seems to run into a problem somewhere below 0.40.
I don't know if it is a genuine problem or a technicality, but you can try this:

res:=dsolve(eval([sys_ode, ics], n = .3787*c),numeric, output=operator,
       continuation = c,maxmesh=2048);
This continuation process starts with c=0 and ends with c=1, i.e. starts with n=0 and ends with n=0.3787 as I have it above. Surely it works by using the previous result as an approximate solution for the present.
Incidentally, you could make that process yourself in a loop by using the optional argument approxsoln=res, where res is the solution found previously. Warning: Don't use output=operator for that process as that currently doesn't work. Just change to output=listprocedure (or the default, procedurelist).

Here is a solution using dsolve/numeric/bvp on the two intervals -0.5..0 and 0..5.
Two parameters are introduced. The parameters are determined by fsolve by requiring smoothness.
You can try the other values in your list a yourself. I show the result for a[1], which is referred to here as A.
The procedure p is primarily for use in finding the parameters, but will also serve to present the result, when output is set at any other name than 'number' (the default).

restart;
sys := diff(f(eta), eta, eta, eta)+f(eta)*(diff(f(eta), eta, eta))+1-(diff(f(eta), eta))^2 = 0, diff(theta(eta), eta, eta)+.72*(diff(theta(eta), eta))+A*(diff(f(eta), eta))+theta(eta) = 0;
bcs := f(-.5) = 0, (D(f))(0) = 1+((D@@2)(f))(0), (D(f))(5) = 0, theta(-.5) = 1+(D(theta))(0), theta(5) = 0;
A:=-.5:
p:=proc(f1,th1,{output::name:='number'}) option remember; local res1,fvals,thvals,res2;
   res1:=dsolve({sys,f(-.5) = 0,D(f)(0) = f1, (D@@2)(f)(0)=f1-1,theta(-.5) = 1+th1,D(theta)(0)=th1},numeric,:-output=listprocedure);
   fvals:=subs(res1,[seq(diff(f(eta),[eta$i]),i=0..2)])(0);
   thvals:=subs(res1,[seq(diff(theta(eta),[eta$i]),i=0..1)])(0);
   ##
   res2:=dsolve({sys,f(0)=fvals[1],D(f)(0) = fvals[2],theta(0)=thvals[1],D(f)(5) = 0,theta(5) = 0},numeric,:-output=listprocedure);
   if output='number' then
      [fvals[3]-subs(res2,diff(f(eta),eta$2))(0),thvals[2]-subs(res2,diff(theta(eta),eta))(0)]
   else
       res1,res2
   end if
end proc;
p1:=proc(f1,th1) p(_passed)[1] end proc;
p2:=proc(f1,th1) p(_passed)[2] end proc;

p(.3,-.2); #Just testing the procedure p
par:=fsolve([p1,p2],[.3,-.2]);
res1,res2:=p(op(par),output=xxx);
plots:-display(plots:-odeplot(res1,[[eta,f(eta)],[eta,theta(eta)]]),plots:-odeplot(res2,[[eta,f(eta)],[eta,theta(eta)]]));
plots:-display(plots:-odeplot(res1,[[eta,diff(f(eta),eta)],[eta,diff(theta(eta),eta)]]),plots:-odeplot(res2,[[eta,diff(f(eta),eta)],[eta,diff(theta(eta),eta)]]));
plots:-display(plots:-odeplot(res1,[[eta,diff(f(eta),eta,eta)]]),plots:-odeplot(res2,[[eta,diff(f(eta),eta,eta)]]));