I shall consider y1 just a third variable until I hear otherwise (see my comment above).

Furthermore I will not be using the name D for the operator, since it is used by Maple already.

Q:=f->unapply(convert(diff(f(x,y,y1),x)+y1*diff(f(x,y,y1),y)+y*y1/x*diff(f(x,y,y1),y1),D),x,y,y1);

Example:

L:=(x,y,y1)->x^2*sin(y)+37*y1;

Q(L)(s,t,u);

Q(L)(8,Pi/2,-6);

## If you prefer you can use a more explicit procedural form, which also has the advantage that x,y,y1 can be declared local to Q:

Q:=proc(f) local expr,x,y,y1;
  expr:=diff(f(x,y,y1),x)+y1*diff(f(x,y,y1),y)+y*y1/x*diff(f(x,y,y1),y1);
  expr:=convert(expr,D);
  unapply(expr,x,y,y1);
end proc;

### Comment:
One could have expected the more direct definition
Q:=f->unapply(D[1](f)(x,y,y1)+y1*D[2](f)(x,y,y1)+y*y1/x*D[3](f)(x,y,y1),x,y,y1);

would work, but it doesn't. Thus the extra step via diff.

You have to realize that (-1)^(1/3) is exp(I*Pi/3)=1/2+1/2*I*sqrt(3), i.e. it is the principal value of the cube root of -1. Thus it is imaginary.
To get the realvalued cube root use surd:

plot(surd(x,3),x=-1..1);
#Or just use convert:
plot(convert(x^(4/3)*sin(1/x),surd),x=-1..1);

See
?surd

This is taking from the example you gave in

http://www.mapleprimes.com/questions/202124-Defining-A-Function-Of-A-PreDefined#comment222568

restart;
n:=7:
g_temp := add((t[i]+1)*ln(t[i]+1)-t[i]*ln(t[i]), i = 1 .. n);
g := unapply(g_temp, seq(t[i], i = 1 .. n));
##First version is cleaner than the second (in my opinion):
p:=proc(x::Vector) local N; global g,step;
  N:=op(1,x);
  x+step*Vector(N,i->D[i](g)(seq(x[j],j=1..N)))
end proc;
p(<1,2,3,4,5,6,7>);
## This uses the arrow notation:
NS:=(x::Vector)->x+step*Vector(op(1,x),i->D[i](g)(seq(x[j],j=1..op(1,x))));
NS(<1,2,3,4,5,6,7>);

You could use this procedure mp instead of modp:

mp:=(x,p)->x-floor((x/p))*p;

odeplot(soln, [mp(theta(t), 2*Pi), omega(t)/`&omega;H`], 0 .. n*Th, labels = ["&theta;(t)","&omega;(t)/&omega;H"], axes = boxed, style = point, size = [.25, .75])

Another version using the builtin trunc:

mpt:=(x,p)->x-`if`(x>=0,trunc(x/p),trunc(x/p)-1)*p;

I vaguely remember a discussion about this on MaplePrimes concerning which version of an extension of modp to reals is better/faster.

Certainly you cannot have an inequality as an initial condition. The error message actually says that, but in a positive way: it tells you what you can have, i.e equations or expressions. Giving an expression 'expr', say, will be interpreted as giving the equation expr=0.

Does this do it:

display(seq(HyperionOrbit(i, `&omega;H`*i, 1, 1), i = -3 .. 3), view = [-3 .. 3, default]);

It is OK in Maple 18.
A look at the code in Maple 18 and Maple 2015 shows that the code has been changed a bit.
In Maple 2015 there is a puzzling line 4:

hor := evalb(`if`(orientation <> (':-default'),orientation = (':-horizontal'),horizontal = true and vertical = false));

which seems to me to return false if orientation=default, because vertical is per default true.

The value of hor is used at the end of the code

  88     if hor then
  89         vp[i] := [-y, x]
           else
  90         vp[i] := [x, y]
           end if;

A workaround is to give the orientation:

DrawNetwork(N,orientation=horizontal);

##Another one (not so nice).
DrawNetwork(N,horizontal=true,vertical=false);
DrawNetwork(N,horizontal,vertical=false); #Shorter, but slightly confusing.

The code would have worked as intended if the default for 'vertical' had been vertical::truefalse := not horizontal

I corrected your equations some and used numerical solution with the parameters option:

restart;
xm:=t*vm; #Position of master
ode1:=diff(f1(t),t)=(xm-f1(t))/sqrt((xm-f1(t))^2+f2(t)^2)*vd;
ode2:=diff(f2(t),t)=-f2(t)/sqrt((xm-f1(t))^2+f2(t)^2)*vd;
res:=dsolve({ode1,ode2,f1(0)=0,f2(0)=h},numeric,parameters=[h,vm,vd]);
res(parameters=[10,1,5]); #Setting the parameters
plots:-odeplot(res,[f1(t),f2(t)],0..2.1);
plots:-odeplot(res,[f1(t),f2(t)],0..2.1,view=[0..2.1,0..10],frames=25); #An animation: Click and run

Note: When the dog joins its master they have the same position, so the square root is zero. This creates the numerical problem at t = 2.0834...

###  The dog follows along after it joins its master (using piecewise):
restart;
xm:=t*vm;
ode1:=diff(f1(t),t)=piecewise(f2(t)>0,(xm-f1(t))/sqrt((xm-f1(t))^2+f2(t)^2)*vd,vm);
ode2:=diff(f2(t),t)=piecewise(f2(t)>0,-f2(t)/sqrt((xm-f1(t))^2+f2(t)^2)*vd,0);
res:=dsolve({ode1,ode2,f1(0)=0,f2(0)=h},numeric,parameters=[h,vm,vd]);
res(parameters=[10,1,5]);
plots:-odeplot(res,[f1(t),f2(t)],0..3,thickness=3);
plots:-odeplot(res,[f1(t),f2(t)],0..5,view=[0..5,0..10],thickness=3,frames=25);

### Although not a dog owner I tried a third version using events instead of piecewise. Also T(t) is given the time of the happy occasion/collision:
restart;
xm:=t*vm;
ode1:=diff(f1(t),t)=(xm-f1(t))/sqrt((xm-f1(t))^2+f2(t)^2)*vd*b(t)+(1-b(t))*vm;
ode2:=diff(f2(t),t)=-f2(t)/sqrt((xm-f1(t))^2+f2(t)^2)*vd*b(t);
res:=dsolve({ode1,ode2,f1(0)=0,f2(0)=h,b(0)=1,T(0)=0},numeric,discrete_variables=[b(t)::boolean,T(t)::float],parameters=[h,vm,vd],events=[[f2(t)=0,[b(t)=0,T(t)=t]]]);
res(parameters=[10,1,5]);
plots:-odeplot(res,[f1(t),f2(t)],0..3,thickness=3);
res(3);
plots:-odeplot(res,[f1(t),f2(t)],0..5,view=[0..5,0..10],thickness=3,frames=25);

Replace the assignment to u(k) with an assignment to some other local variable, e.g. v:

v:=max(exp((0.5*(k_r+1)^2)*0)*(exp(0.5*n*dx*(k_r+1))-exp(0.5*n*dx*(k_r-1))),0);
u:=[op(u),v];

Your u is not an array, it is a list. If you do want an array you could do:

am:=proc(dx,Nminus,Nplus)
local k_r,n,a2,u;
k_r:=1/2;
a2:=a/2;
u:=Array(Nminus..Nplus);
for n from Nminus to Nplus do
u[n]:=max(exp((0.5*(k_r+1)^2)*0)*(exp(0.5*n*dx*(k_r+1))-exp(0.5*n*dx*(k_r-1))),0);
end do;
u
end proc;

Or with indexing from 1:

am:=proc(dx,Nminus,Nplus)
local k_r,n,a2,u;
k_r:=1/2;
a2:=a/2;
u:=Array(1..Nplus-Nminus+1);
for n from Nminus to Nplus do
u[n-Nminus+1]:=max(exp((0.5*(k_r+1)^2)*0)*(exp(0.5*n*dx*(k_r+1))-exp(0.5*n*dx*(k_r-1))),0);
end do;
u
end proc;

You are not using numerical integration. You are computing the symbolic answer (using int) and then you are applying evalf.
Try replacing int by Int. Then you will have numerical integration.
Also try removing evalf, but using int.

Thirdly, try

int((exp(x)*(4420*cos(4)*sin(4)-544*cos(4)^2+148147*exp(-1)-4225*cos(4)-215203)/(71825*exp(1)-71825*exp(-1))-exp(-x)*(4420*cos(4)*sin(4)-544*cos(4)^2+148147*exp(1)-4225*cos(4)-215203)/(71825*exp(1)-71825*exp(-1))+(32/4225)*cos(4*x)^2+(1/71825)*(4225+(2210*x-6630)*sin(4*x))*cos(4*x)+x^2+8434/4225)^2, x = 0 .. 1,method=_RETURNVERBOSE); ## I learned that option from acer.

(The only successful methods FTOC and FTOCMS are more easily seen by wrapping the whole thing in evalf).

The order and/or the arrangement of the terms in the exact result will influence a subsequent use of evalf, i.e. the rounding errors will differ.

In Eqv[i] there is a missing (t) after v[i] in Av*(1-r[i]^2)*v[i].

The parameter omega[st] is given no value.

The initial conditions for the v's are all zero. It follows from the odes for the v's that they will remain zero (by uniqueness).

I tried:
ic:={u[1](0)=0.8, v[1](0)=0.1,u[2](0)=0.8, v[2](0)=0.1,u[3](0)=0.8, v[3](0)=0.1,u[4](0)=0.8, v[4](0)=0.1};
and used omega[st]=50.

sys:=map(op,eval({seq(EqSys[i],i=1..4)},{params}));
res:=dsolve(eval(sys,omega[st]=50) union ic,numeric);
plots:-odeplot(res,[u[1](t),v[1](t)],0..1);

plots:-odeplot(res,[seq([t,v[i](t)+i/10],i=1..4)],0..1,thickness=2);

It seems that this option is not available in odeplot:

restart;
ode:=diff(x(t),t)=x(t);
res:=dsolve({ode,x(0)=1},numeric);
plots:-odeplot(res,[t,x(t)],0..5,sample=[seq(0..5)]);

Error, (in plot/options2d) unexpected option: sample = [0, 1, 2, 3, 4, 5]

Try

HyperionOrbit(.5, 1.8);

You had a name ωH multiplied onto 1.8. That wouldn't make any sense.

K has zeros in its diagonal, thus I don't see any problem with this:

K:=Matrix([<0, -1, 1, -1>,<-1, 0, -1, 1>,<-1, 1, 0,-1>,<1, -1, -1,0>]);
v:=Vector(4,symbol=a);
w:=K.v;

You forgot t in vyH in Eqns. It should be diff(vyH(t), t).
You prevented yourself from discovering that because you put a colon after the assignment to Eqns.