You clearly have too many boundary conditions, so I interpret your problem as one of determining omega so that the extra boundary condition is satisfied.
You will have to use numerical solution.
(Pi should be spelled like that, NOT pi).

restart;
ode:=diff(y(x), x, x) = -(8*omega*(1-exp(-8*x))*exp(-8*x)/(2*x^2)-Pi^2/(32*x))*y(x);
bcs:=y(0)=0,D(y)(0)=1,y(1/2)=0;
dsolve(ode); #We see that no symbolic solution is to be expected
res:=dsolve({ode,bcs},numeric,method=bvp[middefer]);
plots:-odeplot(res,[x,y(x)],0..1/2);
res(0); # You see that omega=0.805476451353708

Three various ways using elementwise operations (~) .

L:=[-3, -2.5, 0, 1, 3/2, Pi, sqrt(5)];
type~(L,nonnegint); #Just the results in order
L=~type~(L,nonnegint); #Equations
zip(`[]`,L,type~(L,nonnegint)); #Lists

restart;
with(LinearAlgebra):
interface(rtablesize=infinity); #Default 10
A:=RandomMatrix(12,datatype=float); #Example
L,V:=Eigenvectors(A):
L; #The eigenvalues
V; # The columns are eigenvectors corresponding to the eigenvalues, in order

Is this what you had in mind?

http://mapleprimes.com/questions/204592-DEplot-Help-For-Plotting-An-ODE-Solution#answer217446

There is another by the same OP:

http://mapleprimes.com/questions/204603-Conditional-ODE-Simplest-Stickslip-Problem

Your problem was probably that you didn't use exp. e^t won't work, it should be exp(t).
Next problem is the comma: it should be period.
Thus the syntax is:

eq:=60=18+69*exp(-0.0491*t);
solve(eq,t);

In view of the definition of a^b for b irrational:

a^b = exp(b*ln(a))

with an appropriate branch of the logarithm used (in Maple the principal branch) we get

evalf(-sqrt(2)..-exp(-1)); #Your range
eval((-2.0)^x, x=-0.987654321);
eval(exp(x*ln(-2.0)),x=-0.987654321); #Same thing

So what did you expect to see plotted?

If s is always a list of one list of equations as is s in your example, then you can just do
J1 := eval(Ja, op(1,s));

or (under the assumption mentioned) you can shorten that to

J1 := eval(Ja, op(s));

Put the files maple.hdb, maple.ind, maple.lib in a folder, called RIFfolder (or whatever) on your computer, say at the location E:/RIFfolder

When you start a worksheet, then start with redefining libname:

libname:="E:/RIFfolder", libname;

After that you should be able to use whatever is in Rif. If it is a package you need to do the ususal
with(Rif);

You need to use a midpoint method as the error message says.
If in addition you raise 'abserr' you get a solution:

dsol5 := dsolve( dsys3, numeric,method=bvp[middefer],abserr=1e-2);
plots:-odeplot(dsol5,[x,w(x)],0..1);

Your first system is a mix of a pde and two odes.
That must be split.
Certainly the method you intend to use 'method = DuFortFrankel' cannot be used for a system even of pdes.

Your second system has highly problematic ibcs.
You write them as:

IBC3 := {diff(Phi(x, (1/2)*h), z)+2*(diff(w(x, (1/2)*h), x, x)) = 0, diff(Phi(x, (-h)*(1/2)), z)+2*(diff(w(x, (-h)*(1/2)), x, x)) = 0, -7*(diff(w(0, z), x, x))+2*(Phi(0, (1/2)*h)-Phi(0, (-h)*(1/2))) = 0, -7*(diff(w(L, z), x, x))+2*(Phi(L, (1/2)*h)-Phi(L, (-h)*(1/2))) = 0, Phi(0, z) = 0, Phi(L, z) = 0, u(0, z) = 0, u(L, z) = 0, w(0, z) = 0, w(L, z) = 0};

When writing ibcs with derivatives you must use the D-notation.
To see one problem with your present formulation just try:

diff(Phi(x, (1/2)*h), z);
#Since there is no z in (Phi(x, (1/2)*h) this will evaluate to 0.
Instead you should write:
D[2](Phi)(x,h/2)

There are several of those.
There may be more problems.

Do you mean the incomplete beta function? The usual one has only 2 arguments.
The incomplete beta function is

B:=int(t^(a-1)*(1-t)^(b-1),t=0..x);

and Maple returns a result.

So if I understand your intention then your beta is
b:=eval(B,{x=sqrt(t),a=1/5,b=1/2});

Next question: By beta^(-1) do you mean the inverse or the reciprocal?

If you mean the reciprocal then:
plot(1/eval(b,t=f/2),f=0..2);

If you mean the inverse then use a parametric plot:
plot([2*b,t,t=0..1]); ##This line is EDITED: See explanation in the comment "Correction" below.

The revised version suffers from the same problem as in

http://www.mapleprimes.com/questions/205212-Problem-With-Dsolve-

Here and there you have to choose between two systems before dsolve (or DEplot) can be expected to work:

restart;

ode1:=a(t)*(diff(a(t), t))+c(t)*(diff(c(t), t))=0; #rhs?
ode2:=a(t)*(diff(a(t), t))+a(t)*(diff(a(t), t))*b(t)*(diff(b(t), t))*c(t)*(diff(c(t), t))=3*t^2; #rhs?
ode3:=a(t)*(diff(a(t), t))+a(t)*(diff(a(t), t))*c(t)*(diff(c(t), t))+a(t)*(diff(a(t), t))*b(t)*(diff(b(t), t))*c(t)*(diff(c(t), t))=2*t^3; #rhs?

solve({ode1,ode2,ode3},diff({a,b,c}(t),t));
sys1,sys2:=allvalues(%); #Two different systems
DEtools[DEplot](sys1,[a(t),b(t),c(t)],t=0.5..1.4,[[a(1)=1,b(1)=2,c(1)=3]],scene=[t,b(t)]);
DEtools[DEplot3d](sys1,[a(t),b(t),c(t)],t=0.5..1.4,[[a(1)=1,b(1)=2,c(1)=3]],scene=[a(t),b(t),c(t)]);

As tomleslie points out you must leave out one of the boundary conditions when k[1] = 0.

The one that you may want to leave out can be found by using a try... catch ... end try construction:

for b in [bcs] do
  BCS:={bcs} minus {b};
  try
    res[b]:=dsolve(eval({ode1,ode2} union BCS, E union {k[1]=0} ) , numeric, method = bvp[midrich], abserr = 1.10^(-8));
  catch:
    print(b,"failed");
  end try
end do;
#You will see that 4 failed, but 2 didn't:
Leaving out either D(f)(0)=1 or D(f)(N)=0 gave results.
plots:-odeplot(res[D(f)(0)=1],[eta,theta(eta)],0..10);
plots:-odeplot(res[D(f)(0)=1],[eta,f(eta)],0..10); #weird result
plots:-odeplot(res[D(f)(N)=0],[eta,theta(eta)],0..10);
plots:-odeplot(res[D(f)(N)=0],[eta,f(eta)],0..10);

# You must decide which you want. I would choose the second.
# INCIDENTALLY: Don't use k both as k[1] and also as a loop variable. Choose another name for the loop variable, like 'i'. Using k and k[1] as you are doing could cause strange errors.

As I see it you already have two solutions, just not the one you think:

eval(res);
RES1:=eval(res[(D@@2)(f)(-1)]);
RES1(0);
plots:-odeplot(RES1, [[x, 1000*f(x)]], 0 .. 1);
RES2:=eval(res[(D@@3)(f)(1)]);
RES2(0); #Same omega as for RES1
plots:-odeplot(RES2, [[x, 1000*f(x)]], 0 .. 1);

You haven't tried pdsolve/numeric, it seems.

To point you in the right direction you can start with this:

pde:=diff(u(x,t),x,x)=diff(u(x,t),t,t)+sqrt(abs(sqrt(2)*diff(u(x,t),x)-u(x,t)));
res:=pdsolve(pde,{u(x,0)=-1,D[2](u)(x,0)=1,u(0,t)=-1,u(1,t)=0},numeric,spacestep=.1);
res:-animate(t=5);

#This doesn't look great, but take a look at the help page too.
You didn't provide any initial or boundary conditions, so I chose some rather arbitrarily.
Also I inserted an 'abs' as you will see.