This works:

R:=rand(0.0..1.0); #Produces floats between 0 and 1
myColors := [ seq( RGB ( R() , R(),R() ) , j = 1 .. 20 ) ] ;
plot([$1..20],x=0..1,color=myColors,thickness=10);

Since you asked for a shooting method I shall give one here.
But I think that you should try the method used by dsolve/numeric/bvp first as you also did and since that was successful (after help from Tom Leslie!) I would recommend staying with that.

Please see my earlier comment about odeplot before trying the following shooting method.
f0, f3, th1, ph1 are parameters to be determined.
We shoot from eta = blt (= 5) for the reason mentioned in my comment.

ics:=D(f)(blt) = 0, (D@@2)(f)(blt) = 0, theta(blt) = 0, phi(blt) = 0,f(blt)=f0,(D@@3)(f)(blt)=f3,D(theta)(blt)=th1,D(phi)(blt)=ph1;
### Using the parameters option:
resS:=dsolve({Eq1, Eq2, Eq3,ics},numeric,parameters=[beta,f0,f3,th1,ph1],output=listprocedure);
### We shall need the following procedures:
F,F1,Th,Ph:=op(subs(resS,[f(eta),diff(f(eta),eta),theta(eta),phi(eta)]));
### For simplicity let us fix beta at beta0=1:
beta0:=1: #Example: First beta value
###We need t0 determine f0,f3,th1,ph1 so that the following procedure p returns [0,0,0,0]:
p:=proc(f0,f3,th1,ph1)  option remember;
  resS(parameters=[beta0,f0,f3,th1,ph1]);
  [F(0),F1(0)-1, Th(0)-1,Ph(0)-1]
end proc;
#Since fsolve handles a list of procedures and not a procedurelist (as p is) we define these:
for i from 1 to 4 do cat(Q,i):=subs(ii=i,proc(f0,f3,th1,ph1) p(_passed)[ii] end proc) end do;
### You notice that p is called from each Q, but option remember makes this cheap.
## We shall now cheat and use values found by dsolve/numeric/bvp (see my comment on odeplot):
res[1](blt);
## Using two significant digits we try:
p(0.57,0.0022,-0.066,-0.039);
## Now use fsolve with start values [0.57,0.0022,-0.066,-0.039]:
fsolve([Q1,Q2,Q3,Q4],[0.57,0.0022,-0.066,-0.039]);
## Result: [.57702521617123716, 0.22545383523878347e-2, -0.66166698346482951e-1, -0.39152140841741792e-1]
## To query parameters do
resS(parameters);
## To set parameters do (example):
resS(parameters=[beta0,.5,0.002,-0.06,-0.04]);
## or you can do the following in any order:
resS(parameters=[beta=beta0, f0=.5, f3=0.002, th1=-0.06, ph1=-0.04]);

If I understand you correctly this should do it:

f:=solve(y=ln(x+2),x);
plot(f,y=-1..2,filled); p:=%:
tr:=plottools:-transform( (x,y)->[y,x]):
tr(p);

IntegrationTools:-GetIntegrand is somewhat more general than the simple and good method shown by Kitonum:

A:=a*Int(v(x),x=0..1)+b;
IntegrationTools:-GetIntegrand(A);

# To see the procedure do:
showstat(IntegrationTools:-GetIntegrand);

I remember seeing this done in a simple way:

http://www.mapleprimes.com/posts/40006-An-Easy-Way-To-Plot-The-Region-Between-Two-Curves

plot(x^2-4,x=0..2, filled, color=blue): plottools:-transform((x,y)->[x,y+4])(%);

There are several other ways. Here is another:

plot([x^2,4],x=0..2,filled,color=[white,green],transparency=0);

To avoid getting RootOf to begin with you can use the option 'explicit':

solve({Q1, Q2, Q3}, {x, y, z}, explicit); # or explicit=true

See, ?solve, details

Another way is to set the environment variable _EnvExplicit to true before solving:

_EnvExplicit:=true;
solve({Q1, Q2, Q3}, {x, y, z});

# So why not set that variable to true by default? One reason is that sometimes output would be extremely large and not very telling:
p:=z^4+a*z^3+b*z^2+c*z+d;
solve(p=0,z);

In recent versions of Maple (2015 and 2016) this example works:

restart;
X:=<x(t),y(t)>;
ode:=diff~(X,t,t)=<-y(t),x(t)>;
dsolve(ode);
dsolve({ode,x(0)=1,y(0)=0,D(x)(0)=0,D(y)(0)=1}); #Mixed input ivp
res1:=subs(%,X);
dsolve({ode,<x(0),y(0)> = <1,0>,<D(x)(0),D(y)(0)> =<0,1>}); #Vector input ivp
res2:=subs(%,X);
res1-res2;

Notice that eq1 involves f only. Solve for the highest derivative, i.e. f '''.

solve(eq1,{diff(f(x),x,x,x)});
EQ1:=op(eval(%,{beta = .1, n = .5 }));

Notice that the denominator on the right is -3.0+.225*f(x)^2.
Thus you should expect a problem if f(x) (or its absolute value) takes on values close to (or at)
sqrt(3/0.225) = 3.6514837 ...
That is highly likely, however. Indeed I tried with my own cruder procedure and found the following graph (using the boundary conditions for f only of course):

My procedure uses an equidistant partition and apparently didn't notice the problem. Using this as an approximate solution in dsolve gave the same error as you report:

Error, (in dsolve/numeric/bvp) Newton iteration is not converging
Suppose for a moment that my "solution" is actually close to some actual solution. Then the denominator and the numerator of rhs(EQ1) would have a common zero at f(x)=3.6514837. Plotting the denominator and numerator with odeplot seems to confirm that.

I shall give some preliminary results to work with.
It doesn't appear to me that continuation in your parameter lambda helps. It appears to me that the problem is the denominator mentioned in my comment above (to Carl Love) together with the boundary condition D(f)(etainf)=0.
Thus I tried the following:

restart;
Digits := 15;
with(plots):
Nr := .1; Nb := .3; Nt := .1; Rb := 0; Lb := 1; Le := 10; Pe := 1; ss := .2; aa := .1; bb := .2; cc := .3; nn := 3/2;
##Notice abs in Eq1:
Eq1 := nn*abs(diff(f(eta), eta))^(nn-1)*diff(f(eta), eta,eta)-(nn+1)/(2*nn+1)*eta*(diff(theta(eta), eta)-Nr*diff(h(eta), eta)-Rb*diff(g(eta), eta)) = 0;
Eq2 := diff(theta(eta),eta,eta)+nn/(2*nn+1)*f(eta)*diff(theta(eta), eta)+Nb*diff(theta(eta), eta)*diff(h(eta), eta)+Nt*(diff(theta(eta), eta)^2) = 0;
Eq3 := diff(h(eta), eta,eta)+nn/(2*nn+1)*Le*f(eta)*diff(h(eta), eta)+Nt/Nb*diff(theta(eta), eta,eta) = 0;
Eq4 := diff(g(eta), eta,eta)+nn/(2*nn+1)*Lb*f(eta)*diff(g(eta), eta)-Pe*( diff(g(eta), eta)*diff(h(eta), eta)+diff(h(eta), eta,eta)*g(eta) ) = 0;
etainf:=10:
##Taking your bcs as stated:
bcs := f(0) = ss/Le*(D(h))(0), theta(0) = lambda+aa*(D(theta))(0), h(0) = lambda+bb*(D(h))(0), g(0) = lambda+cc*(D(g))(0), (D(f))(etainf) = 0, theta(etainf) = 0, h(etainf) = 0, g(etainf) = 0;
## Replacing lambda with 1 as you want the final value of the continuation parameter lambda = 1:
bcs1:=op(eval({bcs},lambda=1));
## Removing the trouble temporarily:
bcs2:=remove(has,{bcs1},D(f)(etainf));
## I had immediate success with D(f)(etainf)=0.07, but not with lower values than 0.07:
res:=dsolve({Eq1, Eq2, Eq3, Eq4} union bcs2 union {D(f)(etainf)=0.07}, numeric, approxsoln=[g(eta)=exp(-eta),h(eta)=exp(-eta),theta(eta)=exp(-eta),f(eta)=3*tanh(eta/3)]);
odeplot(res,[eta,f(eta)],0..etainf);
odeplot(res,[[eta,g(eta)],[eta,h(eta)],[eta,theta(eta)]],0..etainf);
## Using my own solver I got a result for D(f)(etainf)=0.001. Using that result as an approximate solution in dsolve gave success with dsolve as well:



It is clearly not satisfactory that one has to use another solver as a stepping stone, but this may provide an inspiration for your further work with this problem.
## Note: I can get down further, at least to D(f)(etainf)=0.00001, and the graphs look the same, which is comforting.

As Carl is pointing out you need an extra (and nonhomogeneous) boundary condition in order to determine omega. Furthermore you should replace omega^2 by some name, e.g. omega2.
By the way: why do you want to solve this by a shooting method?

Anyway, choosing the first condition that came to my mind g(1)=1  I had no problem finding a solution which turned out to be s(x) = 0 and g(x)<>0:

restart;
dsys3:={-0.326905829596411e-2*g(x)-(diff(g(x), x, x))-(diff(s(x), x))*(diff(s(x), x, x))-(4/3)*omega^2*g(x), -s(x)*omega^2-(-0.573628192993074e-1*sin(0.571756792348295e-1*x)-0.163452914798206e-2*cos(0.571756792348295e-1*x))*(diff(s(x), x))-(1.00327307112014*cos(0.571756792348295e-1*x)-0.285878396174148e-1*sin(0.571756792348295e-1*x)-1)*(diff(s(x), x, x))+0.220893539279189e-4*(diff(s(x), x, x, x, x))-(9/8)*(diff(s(x), x, x))*(diff(s(x), x))^2-(3/4)*(diff(s(x), x, x))*(diff(g(x), x))-(3/4)*(diff(s(x), x))*(diff(g(x), x, x)), (D(g))(1)+(1/2)*(D(s))(1)^2 = 0, g(0) = 0, s(0) = 0, (D(s))(0) = 0, ((D@@2)(s))(1) = 0, ((D@@3)(s))(1) = 0};
dsys3a:=subs(omega^2=omega2,dsys3);
sys,bcs:=selectremove(has,dsys3a,diff);
##
res:=dsolve(sys union bcs union {g(1)=1},numeric);
res(1);
plots:-odeplot(res,[[x,g(x)],[x,s(x)]],0..1,thickness=3);
odetest({s(x)=0},sys); # Indeed.
ode_g:=%[2];
## ode_g is linear and bcs reduces to g(0)=0, D(g)(1)=0. Thus the problem can be solved symbolically in the familiar way. Briefly:
ode:=diff(g(x),x,x)+k^2*g(x) = 0;
dsolve(ode); # We must have only sin(k*x) because of g(0)=0.
cos(k)=0; #This must be satisfied because of D(g)(1)=0.
#Thus:
k=Pi/2+n*Pi;
## And so
0.326905829596411e-2+(4/3)*omega2=(Pi/2+n*Pi)^2;
w2:=solve(%,omega2);
eval(w2,n=0); #Same as the omega2 found above.
#########
## Now to find solutions for which  s(x) <> 0 proceed in principle the same way (forget the last part).
## Example:
res:=dsolve(sys union bcs union {s(1)=.1},numeric);
plots:-odeplot(res,[[x,g(x)],[x,s(x)]],0..1);
res(1);
### As you may know these problems in various forms have come up often in MaplePrimes.

You could define a multiplication operator in a space of functions of two variables as follows.
M sends a function of two variables f into the operator of multiplying by f, so M(f) is the mapping u->f*u:

M:=f->(u->f*u);
M((s1,s2)->s1)(u)(x,t); #Multiplication by the function f given by f(x,t) = x for all x and t.
M((x,t)->x)(u)(x,t); # Same thing
##Your example:
v:=D[1];
w:=M((x,t)->x)@D[2];
((v@w)-(w@v))(u)(x,t);

You can use the parameters option.
Leave out the assignment to n.
Then replace the dsolve statement with:
sol := dsolve({gl1, gl2, gl3, gl4, init1, init2, init3, init4}, numeric,parameters=[n]);
### Example of setting n:
sol(parameters=[4000]);
sol(.5);
odeplot(sol, [x(t), y(t)], t = 0 .. 6.7);
### Now for your actual question. Use this simple procedure:
psol:=proc(n) sol(parameters=[n]); sol end proc;
## Then display the curves for n=1500..9000 (in steps of 1500):
display(seq(odeplot(psol(nn),[x(t), y(t)], t = 0 .. 6.7),nn=1500..9000,1500));

A quote from the help page for Numerical Summation:
   "In the case of infinite sums, Levin's u-transform is used."

You can see the procedure used:

showstat(`evalf/Sum/LevinU_Symb/1_inf`);
showstat(`evalf/Sum/LevinU`);

You have a space after PolynomialInterpolation and before ([Points], y). That space is interpreted by the 2D parser as multiplication, which obviously makes no sense.

The matrix options are unchanged by elementwise operations.
Try:
MatrixOptions(DC);
      shape = [], datatype = complex[8], storage = rectangular, order = Fortran_order
MatrixOptions(DX);
      shape = [], datatype = complex[8], storage = rectangular, order = Fortran_order
# Notice in particular that the datatype is complex[8] in both cases.
## Same result for
MatrixOptions(Re~(DX));
## You can give your matrix explicitly the datatype float:
DX1:=Matrix(abs~(DC),datatype=float);
MatrixOptions(DX1);
       shape = [], datatype = float[8], storage = rectangular, order = Fortran_order
## Finally, by using simplify/zero as in Adri van der Meer's answer, you get:
MatrixOptions(simplify(DX,zero)); #or DC, same answer for the options
       shape = [], datatype = anything, storage = rectangular, order = Fortran_order