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.

You got two problems:

1. The initial conditions contradict the odes.
2. You cannot solve uniquely for the derivatives.

restart;
ics:=a(1)=0,b(1)=0,c(1)=0;
ode:=a(t)*c(t)*(diff(c(t), t))+2*a(t)*(diff(a(t), t))*b(t) = 1, a(t)*(diff(a(t), t))*b(t)*(diff(b(t), t))*c(t)*(diff(c(t), t))+a(t)*(diff(a(t), t))*b(t)*(diff(b(t), t))+a(t)*(diff(a(t), t)) = 1, a(t)*(diff(a(t), t))*c(t)+2*a(t)*b(t)*(diff(b(t), t)) = 1;
eval(convert({ode},D),t=1);
eval(%,{ics}); #Contradiction
##
solve({ode},diff({a,b,c}(t),t));

Which of the roots to choose?