@Carl Love To put the system into a form that allows you to solve for the highest derivatives you can differentiate the first of the equations. This won't introduce any higher derivatives in the system:

ode1:=diff(sys[1],x);
solve({sys[2],ode1}, {diff(w(x),x$4), diff(psi(x),x$3)}); # Possible
# Now add a boundary condition to ensure that the new system is equivalent to the old:
bcs1:=eval(convert(sys[1],D),x=1); # Might as well have used x = 0.
# So new boundary conditions are
BCS:=bcs union {bcs1};
nops(BCS);
## To these must be added the extra condition to determine omega2.
##
By doing that I get 8 successes out of 9 tries without any optional arguments to dsolve expect for the output option, i.e. only
dsolve({sys[2],ode1} union BCS union {b = 10^(-5)}, numeric,output=listprocedure) . 

'undefined' is not the result of an error.
Example:
Heaviside(0);
                                  undefined
But 0/0 produces the error
                            Error, numeric exception: division by zero
i.e. it doesn't result in 'undefined'.

I suppose that you have already looked at the help page:

?undefined

@mwahab I don't know much about this (actually close to nothing), but I tried the following where my result RES differs from the result obtained by map(pdsolve, [res2], parameters = {m}) in the latter missing the case m = 0.

restart;
with(PDEtools);
declare(u(t, x))
pde:=diff(u(t, x), t)-u(t,x)^m*(diff(u(t, x), x))-u(t,x)^m-u(t,x)^m*(diff(u(t, x), x, x))-u(t,x)^m*(diff(u(t, x), x, x, x)) = 0;
res:=DeterminingPDE(pde, u(t, x), integrabilityconditions = false);
nops(res);
res2:=casesplit(res, parameters = {m});
map(pdsolve, [res2], parameters = {m});
n:=nops([res2]);
map2(op,2,[res2]); # The &where's
eqs:=map2(op,1,[res2]); #The equations including conditions on m
sys:=map2(select,has,eqs,{_eta,_xi}); #the pdes
M:=map2(remove,has,eqs,{_eta,_xi}); # the m equations
M1:=subs([]=[m=m,m=m],M); #Slightly artificial
M2:=map2(op,2,M1); #Final version of m equations
RES:=seq( [M2[i],pdsolve(eval(sys[i],M2[i]))],i=1..n); #result
map(pdsolve, [res2], parameters = {m});

See the help page for dsolve:
?dsolve

Include any initial conditions to determine the arbitrary constant(s).

Have a look at the help page
?dsolve,numeric

It has examples of ivp as well as bvp problems.

Note: I must agree with vv that this seems to be a case of using a generic formula.
So I converted this into a comment.
Here is a simplified version, where I set t = 0 and a=b=2 right away.

restart;
G:=n->sum(cos(2*Pi/n*j)^2*sin(2*Pi/n*j)^2,j=0..n-1); #The simplified version
Ga:=n->add(cos(2*Pi/n*j)^2*sin(2*Pi/n*j)^2,j=0..n-1); #Version using add
H:=n->sum(combine(cos(2*Pi/n*j)^2*sin(2*Pi/n*j)^2),j=0..n-1); #Version combining before summing
G(n); #Wrong for some n
H(n); #Seems to be better
combine~([seq(G(n),n=1..5)]);
combine~([seq(Ga(n),n=1..5)]);
[seq(H(n),n=1..5)];

##Test of n/8 as the answer:
simplify(combine~([seq(G(n)-n/8,n=1..25)]));
simplify(combine~([seq(Ga(n)-n/8,n=1..25)]));
simplify([seq(H(n)-n/8,n=1..25)]);

@Carl Love His expected output can be obtained by doing

PDEtools:-declare(phi(X));
Diff(conjugate(phi(X)),x1);
                             

I checked that mgear is still present in Maple 7 by doing
interface(verboseproc=2);
eval(`dsolve/numeric/mgear`);

I looked at the help page for dsolve/numeric/mgear in Maple 7 via What's New. There it says:

"The mgear method is obsolete, and is no longer available in Maple. The rosenbrock and lsode methods are avaialble for the numeric solutions of stiff initial value problems."

@olivertwist My code works in Maple 2016 (but it doesn't use Douglas Meade's Shoot package, so you don't need the with(Shoot) line).

Which Maple version do you have?
The eval[recurse] command is only available in more recent versions. If the result of that command (i.e. sys) is not free of variables with p added to them (as in fp) then you must have an older version.
In older versions you can do the following (which also works in more recent versions):

SBS:=fppp(eta) = diff(fpp(eta), eta), fpp(eta) = diff(fp(eta), eta), fp(eta) = diff(f(eta), eta),gp(eta) = diff(g(eta), eta), mp(eta) = diff(m(eta), eta), np(eta) = diff(n(eta), eta); #The order is important!!!
sys:=subs(SBS,sys1); #The first substitution is applied first, etc.

#####
I'm puzzled by your remark "I solved the problem using desolve without converting the system into first order one (RKF45)".
If you used dsolve with boundary conditions given at two points (here eta=0, eta=blt=1) then you were not using RKF45 as that method is for initial value problems only.

@olivertwist The best way to find out why something is needed is to try with and without.
In this case (for me at least) I had to isolate the derivatives for shoot to work, which meant passing from ODE1 to ODE2.
If I'm right, this simply means that the shoot procedure was written under the assumption that the derivatives are isolated on the left.

@Carl Love Yes, certainly. I just tried your code in Maple 8. It worked.
The copyright is
op(3,eval(`dsolve/numeric/bvp`));
    `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`
## I have also tested the code in Maple 6 and 7. `dsolve/numeric/bvp` doesn't exist in Maple 6, but does in Maple 7, and your code works in Maple 7, but not in Maple 6.

@tomleslie Yes, I missed that one.

@Doug Meade I have been doing what you describe at the end. I use dsolve with the parameters option and fsolve.
For the system in this question your shoot was roughly fifty times faster than mine. I attribute that to my use of fsolve.
Note: I just added LSSolve as a solver. That makes my procedure go about as fast as yours at least in this example.

@vv I found the horses in the subfolder 'images', so I needed:
X:=Read(cat(kernelopts(mapledir),"/data/images/fjords.jpg")):

@Carl Love Yes, I can see that I (among other people) made a comment to that question. I still have a copy.