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?

As tomleslie has pointed out this can be solved by just using dsolve. I think he is right in suspecting that you are not going to like the very complicated answer. The problem is solving a 5th degree polynomial. This is well known to be impossible to do in general (meaning getting a formula for the solution).

I cannot get around this either (of course), but can maybe throw some light on your system.
By making a change of variables r=exp(t) you get a system with constant coefficients.
This system cannot be solved for the highest derivatives (order 2). Thus it is somewhat special.
Anyway here is what I did.
Your 'odesys' is just a copy of what you wrote. No attempt made to make the input look nice.

restart;
odesys := {-(1/2)*(-2*(diff(lambda[1](r), r, r))*phi[0]*r^2+2*(diff(lambda[1](r), r, r))*r^2+(diff(alpha[1](r), r, r))*r^2+2*sigma[1](r)-2+4*(diff(lambda[1](r), r))*r+4*(diff(alpha[1](r), r))*r+2*(diff(sigma[1](r), r))*r)/r^2 = 0, -(1/2)*(4*(diff(lambda[1](r), r, r))*r*phi[0]+6*(diff(alpha[1](r), r, r))*r*phi[0]+8*(diff(lambda[1](r), r))*phi[0]+12*(diff(alpha[1](r), r))*phi[0])/r-(1/2)*(4*(diff(lambda[1](r), r, r))*r*w+12*(diff(lambda[1](r), r))-4*(diff(alpha[1](r), r))+8*(diff(lambda[1](r), r))*w+6*(diff(lambda[1](r), r, r))*r+6*(diff(alpha[1](r), r, r))*r)/r = 0, (1/4)*r*(2*(diff(r, r, r))+2*(diff(sigma[1](r), r, r))*r+2*(diff(alpha[1](r), r))-2*(diff(beta[1](r), r))+4*(diff(sigma[1](r), r)))-(diff(lambda[1](r), r))*r*phi[0]+(1/4)*r*(2*(diff(alpha[1](r), r, r))*r+4*(diff(lambda[1](r), r, r))*r+8*(diff(lambda[1](r), r))+4*(diff(alpha[1](r), r))) = 0, (1/4)*(64*(diff(lambda[1](r), r))+44*(diff(alpha[1](r), r))+8*(diff(r, r, r))*w+16*(diff(lambda[1](r), r))*w+32*(diff(lambda[1](r), r, r))*r+22*(diff(alpha[1](r), r, r))*r)/r+(1/4)*(4*alpha[1](r)-4*beta[1](r)+4*(diff(r^2, r, r))+12*(diff(sigma[1](r), r))*r-4*(diff(beta[1](r), r))*r)/r^2-(alpha[1](r)-sigma[1](r))/(w^2*r^2) = 0};
###Use new names for the dependent variables as functions of the new independent variable t:
fnchange:=remove(has,indets(odesys,function),diff)=~{A(t),B(t),L(t),S(t)};
#The new system
sys:=simplify(PDEtools:-dchange({r=exp(t)} union fnchange,odesys,[t,A,B,L,S]));
#Notice the exp(-2*t) in the first 3 equations. They can be removed:
ode1:=collect(simplify(exp(2*t)*sys[1]),diff);
ode2:=collect(simplify(exp(2*t)*sys[2]),diff);
ode3:=collect(simplify(exp(2*t)*sys[3]),diff,factor);
ode4:=sys[4];
##The system {ode1,ode2,ode3,ode4} has constant coefficients.
solve({ode1,ode2,ode3,ode4},diff({A(t),B(t),L(t),S(t)},t,t)); #Cannot solve for the highest derivatives!
##Try this, but the output is huge, notice the RootOf
dsolve({ode1,ode2,ode3,ode4},{A(t),B(t),L(t),S(t)},method=laplace);
##To gain some insight it helps to look at the result for some concrete values of w and phi[0]:
SYS:=eval({ode1,ode2,ode3,ode4},{w=1,phi[0]=1}); #Special case
var:={A,B,L,S};
inits0:=var(0)=~cat~(var,0);
inits1:=(D~(var))(0) =~cat~(var,1);
###No result comes out of the following, which surely reflects the special nature of the system, i.e. not being possible to solve for the highest derivatives.
#
dsolve(SYS union inits0 union inits1,{A(t),B(t),L(t),S(t)},method=laplace);
#So we try leaving out inits1
res1:=dsolve(SYS union inits0,{A(t),B(t),L(t),S(t)},method=laplace);
## and notice that D(B)(0) is missing, the other derivatives are there.
indets(res1,specfunc(D(B))); #D(B)(0) missing
inits1;
res2:=dsolve(SYS union inits0 union inits1 minus {inits1[2]},{A(t),B(t),L(t),S(t)},method=laplace);
indets(res2,name);
### Is the result useless? Not at all:
RES:=eval(res2,(indets(res2,name) minus {t,_alpha1})=~1); #Setting all initial conditions to 1: just a simple example
ABLS:=subs(RES,[A,B,L,S](t)); #The list of results
eval(ABLS[1],t=1.0); #Check to see if it evaluates fast enough to attempt a plot: yes
plot(ABLS,t=0..5);
plot(eval(ABLS,t=ln(r)),r=0.5..1.5); #Plot as functions of r