restart;
with(plots):
plt := inequal({sqrt((1/2)*x) < y and sqrt(Pi/(2*x)) < y and y < sqrt(x) and y < sqrt(Pi/x)}, x = 1 .. 3, y = 0 .. 2);
display([plt,
  textplot([
    [1.8,0.7, typeset('sqrt(x/2)')],
    [2.5,0.6, typeset('sqrt(Pi/x^2)')]
  ])
]);

In several places in your code you have used square brackets [ ] for grouping.  But square brackaets are not the same as parentheses ( ) in Maple or any other programming language.  Here is the fixed code:

restart;
ddesys := {diff(pred(t), t) = pred(t)*(prey(t)/(prey(t)+10)-2/3), diff(prey(t), t) = prey(t)*(2-(1/20)*prey(t-tau))-prey(t)*(pred(t)/(prey(t)+10))-10, pred(0) = 1, prey(0) = 1};
dsn := dsolve(eval(ddesys, tau = 0), numeric);
plots[odeplot](dsn, [[t, prey(t), color = red], [t, pred(t), color = blue]], t=0 .. 0.3, legend = [prey, pred], labels = [t, ""]);

In the last line I have replaced your plotting range of 0..300 with 0..0.3.  When you execute the code, you will see that the prey population drops down to zero and then becomes negative at time t=0.1.  This indicates a flaw in your model.  See how you can remedy that.

Your solution seems to be correct in general outline although I have not checked the details.  Part of the complication in what you have presented is due to the need to handle the internal forces within the rod.

Lagrangian mechanics (analytical mechanics) introduces a method for deriving a system's equations of motion without the need to deal with the internal forces, and that significantly simplifies the analysis.  Here I illustrate how to solve your problem through the Lagrangian formulation.

Download thrown-bola.mw

There is no expression for a "general solution" to your equation. Normally you would want to solve a PDE along with some boundary conditions which you have not specified.  If you are happy with just any particular solution, consider this one:

mysol := C0 + C1*r*cos(theta) + C2*r*sin(theta);

To verify that this is indeed a solution for any three constants C0, C1, and C2, do

simplify(subs(a(r,theta)=mysol, pde));

where "pde" stands for your PDE.

Producing a good-looking graph with implicitplot3d() requires a fine mesh which is not the most efficient way of plotting a surface and I avoid that whenever possible.  In your case each of the three surfaces may be expressed as a simple function of w, so why not do it with plot3d()?

solve~([f=0, g=0, Op=0], w);
plot3d(%, x=0..2, z=0..2, view=0..2, plotlist, color = ["Red", "Blue", "Green"]);

I have attached an edited version of your worksheet.  It does not solve your problem but contains comments and pointers that can be useful for making progress.

HAM_build-comments.mw

The program still works in Maple 2018.  All I did was to comment out a few lines near the top of the file exam_wkptest.mws which specify a hard-coded path to the library.  Here is the edited version:

restart:
#libname:=`i:\\ptest`,libname;
#march(create,libname[1],100);
#read`i:\\ptest\\wkptest_cpc`;
read "wkptest_cpc";
with(wkptest);

This assumes that the files exam_wkptest.mws and wkptest_cpc are placed in the same folder/directory.

restart;
with(plots):
L := ListTools:-Rotate([$-10..10], 11);
animate(plot,[m*x,x=-10..10],m=L,
    scaling=constrained, view=[-1..1, -1..1]);

You need a couple of fixes to get your code to work.

First, we note that your system of equations is of the 3rd order in u and 4th order in w.  Therefore you need to supply seven boundary conditions but you have supplied only six.  Check what is missing.  Perhaps D(u)(1)=0?

Second, your ode1 is of the second order in u and 3rd order in w.  Differentiate it once in order to put it on a comparable footing with ode2.

Thus, we do:

ode1_new := diff(ode1,r);
ics := ics, D(u)(1)=0;   # change as needed
dsolve({ode1_new,ode2} union {ics},{w(r),u(r)},numeric);

You are attempting to solve the set of three differential equations Eq1, Eq2, Eq3 on the interval [−d, l] but that's not what you really want. What you want is to solve Eq1, Eq2, Eq3 on the three consecutive intervals [−d, 0], [0, p], and [p, l], and connect their solutions through the prescribed interface conditions. In the attached worksheet I show how to do that. Here is the resulting graph of the temperature on the overall interval [−d, l]. The middle interval [0, p] is too tiny to be visible in this graph since p=0.001.

Worksheet:  zz.mw

You have

Phi2:=(t)->Array([seq(...

The error arises when you use Phi2 as if it were a Vector, but Arrays and Vectors are not the same thing in Maple.

Additionally, you should know that if u and v are vectors of length n each, then to calculate their dot product you need to multiply the transpose of u by v, which is written "u^+ . v" (without the quotes).   For that reason,  if A is an nxn matrix, then the bilinear form is "u^+ . (A . v)".  In several places in your code you have written that as the equivalent of u . A . v which is another source of error.

Advice:

Before writing a large and complex code, test small fragments of it elsewhere and only then incorporate that in your final program.  It's easier to catch errors that way.  That's what experienced programmers do.

It's simple.  Just do
 

pde := diff(u(x,t),t,t)-c^2*(diff(u(x,t),x,x)) = 0;
pdsolve(pde);

The result is u(x,t) = _F1(c*t+x) + _F2(c*t-x), where _F1 and _F2 are arbitrary functions.  For instance, you may take _F1 to be the sine function, and _F2 to be the cosine function, but certainly you are not limited to those.

Look at the standard expansion of exp(x) in the power series.  Between the terms insert fracttional powers of x with zero coefficients, as in
1 + 0*x^(1/2) + 1/1!*x^(2/2) + 0*x^(3/2) + 1/2!*x^(4/2) + 0*x^(5/2) + 1/3!*x^(6/2) + O(x^(7/2));

Of course the result is the same as the standard expansion -- the fractional powers have no effect.

What you observe as the parameter "a" changes sign is not a bifurcation.  The system's critical point, that is, (0,0), remains the sole critical point before and after the transition.  What you have here is a loss of stability.  The point (0,0) is a stable equilibrium when a<0 and an unstable equilibrium when a>0.

A bifurcation occurs when an equilibrium splits into multuple equilibria.  A Hopf bifurcation occurs when an equilibrium spawns a periodic orbit.  Admittedly, at a=0 you do get periodic orbits in your case, but that is not what is understood by Hopf bifurcation.  Your system is linear.  A bifurcation is a strictly nonlinear phenomenon.