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.

Leave out rhs as that means "right hand side".

modfit3 := a*x^1.5;
 
f := unapply(modfit3, x);

I agree with tomleslie, you will ultimately have to use numerical solution.
There are a few relevant results that you can get symbolically for this autonomous first order system, however. In particular you can find its equilibria (i.e. stationary solutions) and you have a way of determining their stability.

restart;
sys2 := [diff(x(t), t) = r*x(t)*(K-x(t)-a*y(t))/K, diff(y(t), t) = s*y(t)*(L-y(t)+b*x(t))/L-q*y(t)*E(t), diff(E(t), t) = W*((p-tau)*l*y(t)-c)*E(t)];
params := indets(sys2, name) minus {t}); #Plenty of parameters
#The right hand sides:
f,g,h:=op(subs(x(t)=x,y(t)=y,E(t)=E,rhs~(sys2)));
#The equilibria:
EQ:=solve({f,g,h}=~0,{x,y,E});
pts:=seq(subs(eq,[x,y,E]),eq=EQ); # The equilibria as points
# To determine stability of the pts you need:
J:=unapply(VectorCalculus:-Jacobian([f,g,h],[x,y,E]),x,y,E):
#Here are the corresponding jacobian matrices, whose eigenvalues will determine the stability:
seq(J(op(p)),p=pts);

#Even for the system with E(t)=0 (for all t) is it impossible to get a symbolic solution
sysxy0:=eval(sys2[1..2],E=0);
select(r->r[3]=0,[pts]); #Corresponding equilibria
solxy0:=dsolve(sysxy0); #Keeps on going for longer than my limited patience

Note: These kinds of systems have been studied quite a lot.
A study is made considerably easier if you non-dimensionalize the system. That way your set of 12 parameters can be reduced to 5, which is still a lot, but a huge reduction nevertheless. Non-dimensionalization can be done in many ways. If interested I can give you one choice.

W:=(x,y,yp)->y*yp/x;
W(x,y(x),diff(y(x),x));

You could use lexicographic order defined by f below.

restart;
f:=(z1,z2)->evalb(Re(z1)<Re(z2) or Re(z1)=Re(z2) and Im(z1)<=Im(z2));
C:=RandomTools:-Generate(complex(float(range=-9..9)),makeproc=true); #Generates random complex floats
C(); #Example
N:=50:
L:=[seq(C(),i=1..N)]; #A list of N complex floats
LS:=sort(L,f); #The same sorted lexicograhically
#All the points in the complex plane:
plots:-display(seq(plots:-complexplot(L[1..i],style=point,symbol=solidcircle,symbolsize=15),i=1..N));
#This animation runs through all the points from least to largest (in the lexicographic sense):
plots:-display(seq(plots:-complexplot(LS[1..i],style=point,symbol=solidcircle,symbolsize=15),i=1..N),insequence=true);

The weird constructions I mentioned in my comments above are due to one line:

F := simplify(R+N*U(t), 'size');

This should be corrected to:

F := simplify(R+N*U, 'size');

At the point where the command is located, U has not yet been defined. That is no problem. But when it gets defined, it is not defined as a procedure (i.e. function), thus should not be used as such.

1. Replace omega^2 by a name  e.g.  omega2:

   dsys4 := subs(omega^2 = omega2, dsys3);

2. Choose extra_bcs by these steps:

   sys, bcs := selectremove(has, dsys4, diff);
   extra_bcs := {seq(seq(((D@@i)(g))(a), i = 0 .. 1), a = 0 .. 1), seq(seq((D@@i)(s)(a), i = 0 .. 3), a = 0 .. 1)} minus lhs~(bcs);
  
3. Remove method=bvp[middefer]. No need for it.

4. In dsolve remember to use dsys4 instead of dsys3, and replace omega=34 by omega2=6 (e.g.)

Using Digits:=15 I had success for 2 of the extra conditions:

indices(res);
     [(D(g))(0)], [(D(s))(1)]
They corresponded to omega2 values which can be seen by doing
res[D(g)(0)](0); #omega2 = 1.84809921786846
res[D(s)(1)](0); #omega2= 2.30319572201989

restart;
eq:=exp(-h^2/T)*(Int(exp(-x^2/T)*BesselI(0, h*x/T)*x, x = 0 .. 1))/T;
eq1:=eval(eq,T=1);
#The function BesselI( 0, y) is an even function, thus eq1 is an even function of h.
# Try plotting eq1 as a function of h:
plot(eq1,h=0..10);

This plot suggests that eq1 is decreasing for h>=0. If that is correct then the maximal value of eq1 is
value(eval(eq1,h=0));
     1/2-(1/2)*exp(-1)
i.e. less than 0.99 for sure.

#With some effort one may be able to prove that eq1 is decreasing in h for h>=0.
#The value of eq at h=0 is for general T:
value(eval(eq,h=0));
 # which is  1/2-(1/2)*exp(-1/T), i.e. always less than 1/2.

I don't see any differential equation in your question. Do you want to plot

eq:=q(t)=C*exp(-1/2) + sin(2*t) +cos(4*t) ;

?
I had to correct several errors there, so my interpretation of your intent might be wrong.

If you want to plot this, C needs to have a value, say C=3.
Then you can do

C:=3; #There are better ways than assignment, but use this for now.
plot( rhs(eq),t=0..2*Pi); #rhs = right hand side

f:= x -> (x-a)*(x-b); #Function definition using an arrow ->  . Also notice *

 diff(f(x), x);

I managed without too much effort to get down to (and including) NBT=0.4.
With an additional effort I got down to NBT=0.285.

What I did was to comment out two lines like this:
#NBT:=5:  
#N_bt:=NBT;

Then defined Q,Q1,Q2,Q3,Q4 basically as in my comment above, but with an important difference: NBT is made a formal argument of Q, but not of Q1,Q2,Q3,Q4. Then I ran a loop in NBT backwards from NBT=4.9 to 0.2 in steps of 0.1:

Q:=proc(pp2,fi0,HTCC,ti0,NBT) option remember; local res,F0,F1,F2,F3,a,INT0,INT10,INT15,INT20;
#print(pp2,fi0,HTCC,ti0);
if not type([pp2,fi0,HTCC,ti0],list(numeric)) then return 'procname(_passed)' end if:
res := dsolve({subs(phi0=fi0,p2=pp2,HTC=HTCC,eq2)=0,subs(phi0=fi0,eq1=0),eval(eq3,N_bt=NBT)=0,D(u)(1)=0,u(0)=lambda*D(u)(0),(T)(1)=ti0,phi(0)=fi0,D(T)(1)=0}, numeric,method=bvp[midrich],output=listprocedure);
F0,F1,F2,F3:=op(subs(res,[u(eta),phi(eta),T(eta),diff(T(eta),eta)])):
INT0:=evalf(Int(2*F0(eta)*(1-eta),eta=0..1));
INT10:=evalf(Int(2*F0(eta)*F1(eta)*(1-eta),eta=0..1));
INT15:=evalf(Int((2*F0(eta)*(1-eta)*(F1(eta)*rhop*cp+(1-F1(eta))*rhobf*cbf )),eta=0..1));
INT20:=evalf(Int((2*F0(eta)*F2(eta)*(1-eta)*(F1(eta)*rhop*cp+(1-F1(eta))*rhobf*cbf )),eta=0..1));
a[1]:=INT15-10000*pp2;
a[2]:=INT10/INT0-phi_avg;
a[4]:=F2(0)-1:
a[3]:=INT20/INT15:
[a[1],a[2],a[3],a[4]]
end proc;
Q1:=proc(pp2,fi0,HTCC,ti0) Q(_passed,NBT)[1] end proc;
Q2:=proc(pp2,fi0,HTCC,ti0) Q(_passed,NBT)[2] end proc;
Q3:=proc(pp2,fi0,HTCC,ti0) Q(_passed,NBT)[3] end proc;
Q4:=proc(pp2,fi0,HTCC,ti0) Q(_passed,NBT)[4] end proc;

inipt:=[80,0.03,2,2]:
for NBT from 4.9 to 0.2 by -0.1 do
 ResOpt[NBT]:=Optimization:-LSSolve([Q1,Q2,Q3,Q4],initialpoint=inipt);
 inipt:=convert(ResOpt[NBT][2],list);
end do;

sort([indices(ResOpt,nolist)]);
Successful all the way down to NBT=0.4. Fails at NBT=0.3.
###################
Starting another loop with increments of -0.01 instead of -0.1.
inipt:=convert(ResOpt[0.4][2],list);
for NBT from 0.395 to 0.2 by -0.01 do
 ResOpt[NBT]:=Optimization:-LSSolve([Q1,Q2,Q3,Q4],initialpoint=inipt);
 inipt:=convert(ResOpt[NBT][2],list);
end do;
#Successful all the way to NBT=0.285. Failure at NBT=0.275.
#I didn't have success with starting at NBT=0.39 for some reason, so chose NBT=0.395.
#You could attempt a third loop starting with NBT=0.2845 and in increments of -0.005 or whatever happens to work and using initially before the new loop:
inipt:=convert(ResOpt[0.285][2],list);
#which I found to be:
inipt := [68.1893210471432, 0.224788905364089e-3, 3.97466998520601, -.772490643729996]

With

p := display({p1, p2, p3});

you can use

plottools:-getdata(p);