Doesn't the way I showed you previously work?

http://www.mapleprimes.com/questions/201155-Programme-For-Nonlinear-Ode

The error message tells you what to try: a midpoint method.

Since L only has 3 elements you obviously cannot let k go to 10. I changed it to 3.
Also you were overwriting the previous results so I changed to sol[k].
Furthermore I use Array instead of array.

for k to 3 do
  sol[k] := dsolve(eval({bcs1, bcs2, eq1, eq2}, k_1 = L[k]), [f(eta), theta(eta)], numeric,
  output = Array([0., .1, .2, .5, 1.0, 2.0]), method = bvp[midrich], maxmesh = 1024)
end do;

sol[1];
sol[2];
sol[3];
col:=[red,blue,green]:
display(seq(odeplot(sol[i],[eta,f(eta)],0..2,color=col[i]),i=1..3));

The initial conditions ought to be

DEplot({eq1, eq2}, [v1(t), v2(t)], 0 .. 3*Pi, {[init1,init2]}, stepsize = .1);
or
DEplot({eq1, eq2}, [v1(t), v2(t)], 0 .. 3*Pi, {[0, 2, 0]}, stepsize = .1);

However, it is interesting that initial conditions v1(0)=0, v2(0)=0 creates an error. After all, the zero solution v1=0, v2=0 is obvious.
odetest([v1(t)=0,v2(t)=0],{eq1, eq2});

Using no numerics gives the (unique) solution to v1(0)=0, v2(0)=0

dsolve({eq1,eq2,v1(0)=0,v2(0)=0});

Numerical solution directly (not via DEplot)
res:=dsolve({eq1,eq2,v1(0)=0,v2(0)=0},numeric);
res(1.234567);

plots:-odeplot(res,[[t,v1(t)],[t,v2(t)]],0..3*Pi,thickness=5);

Thus there seems to be a DEplot problem.
My guess is that it is not handling starting at an equilibrium point too well.
Notice that the error shows up in `DEtools/DEplot/direction`, thus likely has to do with drawing the arrows. But there will be no direction to point at at the point (v1,v2)=(0,0):
This is confirmed by:
DEplot({eq1, eq2}, [v1(t), v2(t)], t=0 .. 3*Pi, {[v1(0)=0, v2(0)=0]},arrows=none);
and even more so by
DEplot({eq1, eq2}, [v1(t), v2(t)], t=0 .. 3*Pi, {[v1(0)=0, v2(0)=0]},v1=-1..1,v2=-1..1);

In many (most) situations dsolve can find the dependent variable(s) itself. However, in some cases, as in this, it needs to be told:

dsolve({Eq1,bcs},psi(y));

Notice that dsolve/series and dsolve/laplace always need the dependent variable.

Initially two remarks:
1. You have an ode, not a pde.
2. Don't use GAMMA as that is the name of the gamma function in Maple. You can use Gamma instead.

I'm assuming that you are trying to diagonalize the system {deq1,deq2}.
For completeness I included all the lines.
restart;
with(PDEtools): with(LinearAlgebra):
ode := diff(u(t), t, t)+2*Gamma*(diff(u(t), t))+omega^2*u(t) = 0;
deq1 := diff(u(t), t) = v(t);
deq2 := subs(deq1, ode);
dsolve({deq1, deq2}, {u(t), v(t)});
eqns := [rhs(deq1) = lhs(deq1),rhs(deq2) = lhs(deq2)];
y := [u, v]; #
A, b := GenerateMatrix(eqns, y(t));
gnat := Eigenvectors(A);
lambda := gnat[1]; Lambda := gnat[2];
Y := Vector(y);
GenerateEquations(Lambda, [x[1], x[2]], Y);
tr1 := solve(%, {u, v}); #
tr:=subs(x[1]=x[1](t),x[2]=x[2](t),u=u(t),v=v(t),tr1); #Take notice
dchange(tr, {deq1, deq2}, [x[1], x[2]]); #
solve(%,{diff(x[1](t),t),diff(x[2](t),t)}); #

ode1:=diff(y(x),x)=x/2/y(x); #ode satisfied by x^2-2*y^2 = C
ode2:=diff(y(x),x)=-2*y(x)/x; #ode satisfied by the orthogonal trajectories: the product of the right hand sides should be -1.
res1:=dsolve(ode1,implicit);
res2:=dsolve(ode2);
#Thus using either implicitplot or contourplot as proposed by Carl and Kitonum you get:
p1:=plots:-contourplot(x^2-2*y^2, x= -5..5, y= -5..5, contours= [seq(C, C= -5..5)],color=blue):
p2:=plots:-contourplot(x^2*y, x= -5..5, y= -5..5, contours= [seq(C, C= -5..5)],color=red):
plots:-display(p1,p2);
p1:=plots:-implicitplot([seq(x^2-2*y^2=C, C=-5..5)], x= -5..5, y= -5..5, gridrefine=2,color=blue):
p2:=plots:-implicitplot([seq(x^2*y=C, C=-5..5)], x= -5..5, y= -5..5, gridrefine=2,color=red):
plots:-display(p1,p2);

I'm pretty convinced that in fact your problem has no solution.
The short reason is that in order to avoid an actual singularity where 1-xi*f(eta)^2 = 0 you must impose a constraint on the derivatives at such a point.
This means that for your third order eq1 you can only impose 2 boundary conditions, not 3.
Thus if you impose f(0) = 0 and D(f)(0) = 0 then no more conditions are allowed.

I went through the method described for your eq1, but it is lengthy, and since the result is negative, I will present a much simpler example of the same thing here.

restart;
ode:=t^2*diff(x(t),t,t)+t*diff(x(t),t)+t^2*x(t)=0; #Bessel equation of order 0.
eq:=expand(isolate(ode,diff(x(t),t,t))); #Singularity at t = 0.
dsolve(ode);
dsolve({ode,x(0)=1}); #Notice: no arbitrary constant!
plot([BesselJ(0,t),BesselY(0,t)],t=-3..3);
Now a numerical approach that takes into account the constraints posed by the singularity at t=0.
First of all it follows immediately from eq that x ' (0) must be zero for a solution to exist on an interval containing 0 as an interior point.
We expand x(t) in a taylor series  to order 2:
T2:=eval(mtaylor(x(t),t=0,3),{x(0)=x0,D(x)(0)=0,(D@D)(x)(0)=x2});
Insert that into eq:
eval(eq,x(t)=T2);
isolate(%,x2);
Thus we have the following constraint between x0 = x(0) and x2 = x '' (0):
limit(%,t=0);
In this case we can check this with the exact solution found above:
taylor(BesselJ(0,t),t=0);
Since here x(0) = 1, x ' (0) = 0, x '' (0) = -1/2 we have agreement.

Now we modify our equation eq to take the constraint into account.
The right hand side is going to be:
X2:=unapply(piecewise(t=0,-1/2*x0,rhs(eq)),t);
thus the equation is
eqa:=diff(x(t),t,t)=X2(t);
Solving numerivally with x0 as the only parameter that can be varied:
resnum:=dsolve({eqa,x(0)=x0,D(x)(0)=0},numeric,parameters=[x0],output=listprocedure);
Trying it out:
resnum(parameters=[1]);
plots:-odeplot(resnum,[t,x(t)],-3..3);
Comparing to the exact result (which also had x0 = 1) we see that there is excellent agreement:
plots:-odeplot(resnum,[t,x(t)-BesselJ(0,t)],-3..3);
However, the basic point is that we only have one parameter to vary. Thus in general we can impose only one requirement of the type x(t1) = a0 or D(x)(t1) =a1.  
We take an axample where we impose x(-3) = -1.
X,X1:=op(subs(resnum,[x(t),diff(x(t),t)]));
Q:=proc(x0) resnum(parameters=[x0]); X(-3)+1 end proc;
r:=fsolve(Q,2);
resnum(parameters=[r]);
X(-3);
plots:-odeplot(resnum,[t,x(t)],-3..3);

The dependent variable(s) in dsolve should be exactly that: variables.

Try also
dsolve({Eq1,cs});
dsolve({Eq1,cs},f(x));
f(x);

This is the exact same question (with xi instead of pi) as posed on January 29 by kyanoosh.

http://www.mapleprimes.com/questions/200820-How-Do-I-Solve-A-Nonlinear-Differential

This problem appears to be difficult.
No reason to look at eq2. It is eq1 that is the problem.

eq1:=diff(f(eta),eta$3)+(1)/(2)*f(eta)*diff(f(eta),eta$2)-xi*(2*f(eta)*(diff(f(eta),eta))*
(diff(f(eta),eta,eta))+f(eta)^2*(diff(f(eta),eta,eta,eta))+eta*(diff(f(eta),eta))^2*(diff(f(eta),eta$2)))-K*
(diff(f(eta),eta)-1)=0;
#Isolating f ''' :
Eq:=isolate(eq1,diff(f(eta),eta$3));

Now suppose that there is indeed a solution defined for all eta >= 0. Since that solution satisfies f ' (eta) -> 1 as eta - > infinity and f(0)  = 0 certainly f(eta) will attain the value 1/sqrt(xi) for at least one value of eta.
That means that the denominator 1-xi*f(eta)^2 in rhs(Eq) will become zero. Thus you have a singularity to deal with.
That is going to be your worst problem. You need the numerator to be zero when the denominator is, or more precisely, you need the limit rhs(Eq) to exist at any zero of the denominator.
Good luck!

Well your binfinitive (I suppose it is a replacement for infinity, I shall call it inf for short) is really just another parameter as can be seen from a change of variables.
Thus it has no special standing as such. It is just as difficult to deal with as any other parameter might be:

restart;
MUR:=(1-phi)^2.5:
RhoUR:=(1-phi+phi*rho[p]/rho[f]):
RhoCPR:=(1-phi+phi*rhocp[p]/rhocp[f]):
BetaUR:=(phi*rho[p]*beta[p]+(1-phi)*rho[f]*beta[f])/(RhoUR*rho[f])/beta[f]:
dqu3:=diff(h(x),x$1)-RhoUR*BetaUR*T(x);
dqu2:=5*diff(T(x),x$2)+k[f]/k[nf]*Pr*RhoCPR*f(x)*diff(T(x),x$1);
dqu1:=5/(MUR)*diff(f(x),x$3)
+ 2*(diff(h(x),x$1)*x-h(x))
+RhoUR*(3*f(x)*diff(f(x),x$2)-diff(f(x),x$1)^2);
PDEtools:-dchange({f(x)=fy(y),T(x)=Ty(y),h(x)=hy(y),x=y*inf},[dqu1=0,dqu2=0,dqu3=0],[fy,Ty,hy,y]);
sysy:=solve(%,{diff(fy(y),y,y,y),diff(Ty(y),y,y),diff(hy(y),y)});
#The boundary conditions:
bcsy:={Ty(0)=1,Ty(1)=0,fy(0)=0,D(fy)(0)=inf*lambda,D(fy)(1)=0,hy(1)=0};
indets(sysy,name);

If you are talking about ordinary differential equations:

?dsolve,bvp

The main problem is the definition of f. There is a global p hidden in delta. Thus you need unapply:

f:=unapply(BesselJ(0,R*sqrt(-delta*p))*k*p-(R*sqrt(-delta*p))*BesselJ(1,R*sqrt(-delta*p))*(d*p/R + 2*beta*k/(R^2)),p);

By e^x presumably you mean exp(x).
I used E(t) instead of your more elaborate expression.
restart;
eq:=.5*(diff(E(t), t)) = -(-.132+.1*exp(.6*E(t)))*E(t)+0.5e-1;
dsolve(eq); #implicitly given general solution
dsolve({eq,E(0)=1}); #Attempt at solving initial value problem
res:=dsolve({eq,E(0)=1},numeric); #Numerical solution
plots:-odeplot(res,[t,E(t)],0..50);

Try this method first. After that you get off easier by the method at the bottom.
1). Solve the Eq1 bvp problem first.
2). Notice that if f solves that problem then g(x) = -f(-x) (x in -1..0) will solve Eq1 with g(-1)=0, D(g)(-1)=1, g(0)=0,(D@D)(g)(0)=0, D(g)(0)=D(f)(0), and (D@@3)(g)(0) = (D@@3)(f)(0).
3). It follows that we can extend the f found in 1) smoothly by f(x) = -f(-x) for x in -1..0.
4). Now solve the bvp problem for Eq2 with the known option.

Eq1:=diff(f(x),x$4)+(f(x)*diff(f(x),x$3)-diff(f(x),x$1)*diff(f(x),x$2))=0;
Eq2:=diff(f2(x),x$2)+f(x)*diff(f2(x),x$1)=0;
bcsf:=f(0)=0,(D@@2)(f)(0)=0,f(1)=0,D(f)(1)=1;
resf:=dsolve({Eq1,bcsf},numeric,output = listprocedure);
F,F1,F2,F3:=op(subs(resf,[seq(diff(f(x),[x$i]),i=0..3)]));
##Illustrations of smoothness:
plot(-F(-x),x=-1..0):
plot(F,0..1):
plots:-display(%%,%);
plot(F1(-x),x=-1..0):
plot(F1,0..1):
plots:-display(%%,%);
plot(-F2(-x),x=-1..0):
plot(F2,0..1):
plots:-display(%%,%);
plot(F3(-x),x=-1..0):
plot(F3,0..1):
plots:-display(%%,%);
##The extended f:
FF:=x->piecewise(x<0,-F(-x),F(x));
resf2:=dsolve({subs(f=FF,Eq2),f2(-1)=0.5,f2(1)=1},numeric,output = listprocedure,known=FF);
plots:-odeplot(resf2,[x,f2(x)],-1..1);
##############
##Now that we know that f is going to be an odd function we can solve the problem much easier:
bcsfB:=f(-1)=0,(D)(f)(-1)=1,f(1)=0,D(f)(1)=1,f2(-1)=0.5,f2(1)=1;
resfB:=dsolve({Eq1,Eq2,bcsfB},numeric,output = listprocedure);
plots:-odeplot(resfB,[[x,f(x)],[x,f2(x)]],-1..1);

Maybe I misunderstand your problem because it seems obvious to me: You need
alpha>=-(1/2)*S^2/(K+2)
in order to avoid imaginary numbers.
Plotting should not in principle be a problem either since plot just ignores imaginary numbers:
plot3d(eval(lambda1,K=0.12345),alpha=-15..15,S=-15..15);
plot(eval(lambda1,{K=0.12345,S=2.678}),alpha=-15..15);