I shall assume that all variables are positive. You can do with less, but try:

Ex3 := Pi*sqrt((4/mu+4/r)*exp(-r/mu));
simplify(Ex3) assuming positive;
combine(%) assuming positive;

You have eq3= ... and eq4= ... .
You need assignment, so use eq3:= ... and eq4:= ....

I tried dsolve/numeric:

res:=dsolve(ODE union BC1 union {f(0) = 0, fp(0) = .2*fpp(0),h(0) = 0,i(0) = 0,g(0) = 1+.2*gp(0)},numeric);
Why do you want to use a shooting method?
Anyway, if you do, try this:

1. replace ODE by ODE1 given by
ODE1:=solve(ODE,indets(ODE,specfunc(diff)));

2. Replace IC by IC2 given by the explicit initial conditions:
IC2 := {f(0) = 0, fp(0) = gamma1*delta, fpp(0)=delta, g(0) = 1+gamma1*beta, gp(0) = beta, h(0) = 0, hp(0) = beta1, i(0) = 0, ip(0) = beta2, fppp(0) = alpha};

3. Then do:

S := shoot(ODE1, IC2, BC1, FNS, [alpha = .1, beta = .2, beta1 = .3, beta2 = .4,delta=0]);

I tried to insert a print statement in both obj1 and obj2, here only shown in obj1:
 

obj1:=proc(u1,u2)
local z1,s1,s2,z2,res;
global sol1;
sol1('parameters'=[1,0,u1]);
z1 := sol1(.5);
s1:=subs(z1,y1(t));
s2:=subs(z1,y2(t));
sol1('parameters'=[s1,s2,u2]):
z2:=sol1(0.5);
res:=-subs(z2,y2(t));
print(res);
res
end proc;

Then I tried your worksheet. That showed that fdiff in computing Hess1 uses only 18 digits (i.e. the global setting of Digits + 3), whereas in Hess2 twice as many digits are used.Then I tried a third version obj3, where sol1 is made local. That mostly defeats the purpose of using the parameters option in dsolve, I suppose, but here it is:

obj3:=proc(u1,u2)
local z1,s1,s2,z2,res,sol1;
sol1 := dsolve({sys, y1(0) = alpha, y2(0) = beta}, type = numeric, 
'parameters' = [alpha, beta, u],maxfun=0,range=0..0.5):
sol1('parameters'=[1,0,u1]);
z1 := sol1(.5);
s1:=subs(z1,y1(t));
s2:=subs(z1,y2(t));
sol1('parameters'=[s1,s2,u2]):
z2:=sol1(0.5);
res:=-subs(z2,y2(t));
print(res);
res
end proc;

I then got (almost) exactly Hess2, when doing:

Hess3:=Matrix(2,2):
for i from 1 to 2 do 
  for j from 1 to 2 do 
    Hess3[i,j]:=fdiff(obj3,[i,j],u0,workprec=1.0)
  od
od:
Hess3;

## I tried Digits:=10, then I got Hess1 roughly equal to Hess2.
 

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});