Change the dsolve line to:

dsol[i] := CodeTools:-Usage(dsolve(dsys, numeric, continuation = lambda));

You mentioned method RK45 (rkf45?) in your title. Your problem, however, is a boundary value problem and is handled by `dsolve/numeric/bvp`. Method rkf45 is used for initial value problems only.

You can solve the problem by the textbook method of Fourier series.
The first step is a well known trick which is used to reduce the problem to a homogeneous one.

restart;
pde:=diff(u(x, t), t, t)-c^2*(diff(u(x, t), x, x))-cos(x) = 0;
IBC :={u(0, t) = A, u(1, t) = B, u(x, 0) = 0, (D[2](u))(x, 0) = sin(x)};
eval(pde,u(x,t)=a*x+b+cos(x)/c^2);
sol1:=u(x,t)=a*x+b+cos(x)/c^2;
eval(rhs(sol1),x=0)=A;
eval(rhs(sol1),x=1)=B;
solve({%%,%},{a,b});
sol2:=eval(sol1,%);
pde2:=simplify(eval(pde,u(x,t)=v(x,t)+rhs(sol2)));
pdsolve(pde2,HINT=X(x)*T(t));
dsolve(diff(X(x), x, x) = -k^2*X(x));
solve(sin(k)=0,k,allsolutions);
#k = n*Pi
Tsol:=dsolve(diff(T(t),t,t) = -c^2*n^2*Pi^2*T(t));
S:=Sum(sin(n*Pi*x)*(a(n)*cos(Pi*n*c*t)+b(n)*sin(Pi*n*c*t)),n=1..infinity); #Solution for v(x,t)
eval(S,t=0)=-rhs(sol2); #u(x,0)=0
res_a:=a(n)=2*int(rhs(%)*sin(n*Pi*x),x=0..1) assuming n::posint;
diff(S,t); #We can ignore the time independent sol1
eval(%,t=0)=sin(x);
res_b:=b(n)=1/Pi/n/c*2*int(sin(x)*sin(n*Pi*x),x=0..1) assuming n::posint;
S1:=eval(S,{res_a,res_b});
res:=u(x,t)=rhs(sol2)+S1;
plots:-animate(plot,[eval(rhs(res),{A=5,B=9,c=1,infinity=50}),x=0..1],t=0..10,frames=100);

As I said in my comment, you need to rewrite the boundary conditions to get rid of f'''(0) or in the first order systems version v'(0).
Here I do that for the first order system. I do not use a shooting method.
#####
restart;
varsigma := 10;
a := 1;
b := -1;
ss := 2;
Pr := 1;
lambda := 0;
sigma := -2;
ODE := {(diff(t(eta), eta))/Pr+f(eta)*t(eta)-u(eta)*s(eta) = 0, diff(v(eta), eta)+f(eta)*v(eta)-u(eta)^2+lambda*s(eta) = 0, diff(f(eta), eta) = u(eta), diff(s(eta), eta) = t(eta), diff(u(eta), eta) = v(eta)};
## Your boundary conditions:
bcs:={u(varsigma) = 0, s(varsigma) = 0,f(0) = ss, u(0) = sigma+a*v(0)+b*(D(v))(0), s(0) = 1};
## We need to get rid of D(v)(0).
## So we use the ode for v:
ode_v:=diff(v(eta), eta)+f(eta)*v(eta)-u(eta)^2+lambda*s(eta) = 0;
eval(convert(ode_v,D),eta=0);
dv0:=solve(%,{D(v)(0)});
## The rewritten boundary conditions:
bcs2:=eval(bcs,dv0);
res:=dsolve(ODE union bcs2,numeric);
plots:-odeplot(res,[[eta,f(eta)],[eta,s(eta)]],0..varsigma);

The first order systems approach gave answers right away. The original version you posted gave some problems, I guess because of the one nonlinear boundary condition. This condition is also present in the first order version, but just happened not to cause a problem.
I tried using the result from the first order system as an approximate solution for the original. It was comforting that that worked fine:
ode1:=diff(f(eta),eta$3)+f(eta)*diff(f(eta),eta,eta)-diff(f(eta),eta)^2+lambda*theta(eta)=0;
ode2:=1/Pr*diff(theta(eta),eta,eta)+f(eta)*diff(theta(eta),eta)-diff(f(eta),eta)*theta(eta)=0;
## We need to rewrite this:
bcsE:={f(0)=ss,D(f)(0)=sigma+a*(D@D)(f)(0)+b*(D@@3)(f)(0),theta(0)=1,D(f)(varsigma)=0,theta(varsigma)=0};
eval(convert(ode1,D),eta=0);
f30:=solve(%,{(D@@3)(f)(0)});
bcsE2:=eval(bcsE,f30); #The rewritten boundary conditions
## Now we use the result from the first order system:
resM:=dsolve(ODE union bcs2,numeric,output=mesh); #output=mesh is easiest to deal with right now.
## We need to change names and such:
resM2:=subs(s(eta)=theta(eta),t(eta)=diff(theta(eta),eta),u(eta)=diff(f(eta),eta),v(eta)=diff(f(eta),eta,eta),resM);
## Now using resM2 as an approximate solution:
sol:=dsolve({ode1,ode2} union bcsE2,numeric,approxsoln=resM2);
plots:-odeplot(sol,[[eta,f(eta)],[eta,theta(eta)]],0..varsigma); # Same as shown above
##Try first without an approximate solution. You will see problems, which I'm sure is because of the one nonlinear boundary condition.
##############
## Exploring the effect of rewriting the nonlinear boundary condition in the latter formulation:
eqE:=select(has,bcsE2,D(f)(0));
#Choosing to isolate f '' (0):
d20a:=eval(solve(eqE,{(D@D)(f)(0)}),f(0)=ss); # Also explicitly putting f(0)=ss.
d20b:=solve(eqE,{(D@D)(f)(0)}); # Leaving to dsolve to use f(0)=ss.
bcsE3a:=remove(has,bcsE2,D(f)(0)) union d20a;
bcsE3b:=remove(has,bcsE2,D(f)(0)) union d20b;
## Now solving the system with these two formulations of the boundary conditions. We use the same approximate solution in the two cases.
## Using bcsE3a:
resE3a:=dsolve({ode1,ode2} union bcsE3a,numeric,approxsoln=[f(eta)=ss*exp(-2*eta)+ss,theta(eta)=exp(-2*eta)],initmesh=2048);
plots:-odeplot(resE3a,[[eta,f(eta)],[eta,theta(eta)]],0..varsigma); pa:=%:
## Using bcsE3b:
resE3b:=dsolve({ode1,ode2} union bcsE3b,numeric,approxsoln=[f(eta)=ss*exp(-2*eta)+ss,theta(eta)=exp(-2*eta)],initmesh=2048);
plots:-odeplot(resE3b,[[eta,f(eta)],[eta,theta(eta)]],0..varsigma); pb:=%:
plots:-display(pa,pb);  #Surprise! Two solutions.

#Finally (I hope) I tried my own shooting procedure and was without problem able to reproduce the two different solutions.

You don't tell us what f__r and f__h are, but I made up a simple example, where I replaced your large number 9999999 by 'undefined':

#First using 9999999:
restart;
q:= proc (x,y,z)
  if x>0 then x^2+y^2+z^2-1 else 9999999 end if
end proc;
plots:-implicitplot3d('q(x,y,z)',x=-2..2,y=-2..2,z=-2..2);
#Now using 'undefined':
restart;
q:= proc (x,y,z)
  if x>0 then x^2+y^2+z^2-1 else undefined end if
end proc;
plots:-implicitplot3d('q(x,y,z)',x=-2..2,y=-2..2,z=-2..2);

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