implicitplot3d needs a lot of time and memory to produce smooth-looking results, because it's working on a 3-d grid of points.  But if you're using style=contour, most of that is wasted, because you're just plotting a few curves at constant values of z.  The parametric form using plot3d is much more efficient since you're using just a 2-d grid, and so you can get much smoother curves.  Even better might be to plot the contour lines directly using spacecurve.

It's also more efficient to use intersectplot where at least one of the surfaces is given parametrically rather than implicitly.  Thus:

> with(plots):
   display([ seq(spacecurve([x,1/x,z], x=1/3 .. 3, colour=black,thickness=2), z = [$0..5]/2),
                   seq(spacecurve([x, 2*x/(1+2*z),z], x=0 .. 3, colour=black,thickness=2), z = [$0 .. 5]/2),
                   intersectplot(surface([x,2*x/(1+2*z),z],x=0..3,z=0..3),
          surface(y = 1/x, x = 0 .. 3, y = 0 .. 4, z = 0 .. 2.5), color = red, thickness = 2)], axes=normal);

The eggs in the link you provide are not in that form.  The equation y = f(x) would have to be defined piecewise, since these eggs are formed from circular arcs.   In fact the "golden egg" is of the following form:

> yegg := piecewise(x < 0, sqrt(1-x^2), x < sqrt(2), sqrt(4-x^2)-1, sqrt(5-4*sqrt(2)-x^2+2*x));
   plot([yegg, -yegg], x=-1 .. 3 - sqrt(2), colour=red, scaling=constrained);

I don't know what's particularly "golden" about it.

You don't need the a < 2*Pi, but you do need a > 0.  In fact, for a < 0 the result should be -4*a^5*Pi/15.

Do you mean this?

> diff(omega_3(t), t) := alpha_3(t);

Or do you want the analogous result for all omega_n?

select(i -> (M[i] - M[i-1] = H), [$2..nops(M)]);

> plot3d([[x,1/x,z],[x,2*x/z,z]], x=0..4, z=0..4, view=[0..4,0..4,0..4], axes=box, 
    labels=[x,y,z], orientation=[-80,20]);

It's just how floating-point arithmetic works.  The lowest-order digits of x^3 and 2*x, which are the only ones that differ when x is -0.4533976515, -0.4533976516 or -0.4533976517, are lost when you add them, using the default Digits = 10. 

With  x =  -.4533976515,

x^3 = -.09320469696
2*x = -.9067953030
x^3 + 2*x = -1.000000000
which is indeed correct to 10 digits.

With x = -.4533976517

x^3 = -.09320469708
2*x = -.9067953034
x^3 + 2*x = -1.000000000

which is also correct to 10 digits.

Here are the solutions that are functions of x alone.

> with(VectorCalculus); SetCoordinates(cartesian[x,y,z]);
  F:= VectorField(<f1(x),f2(x),f3(x)>);
  eq:= Laplacian(F) + 2*b*c*Curl(F) + c*F;
  dsolve(convert(eq,set));

As your three parameters vary, your vector will (in general) trace out a solid region in space.  Maple doesn't have facilities for graphing this directly.  You can plot points (with pointplot3d), curves (spacecurve) or surfaces (plot3d).
For example, you might use plot3d for the boundaries of your region (which might correspond to theta1 = -Pi or Pi, theta2 = -Pi or Pi, theta3 = -Pi or Pi).

I tried it with pdsolve, but Maple couldn't solve it ina reasonable time (one try got a memory allocation error, another I stopped with Time 4954 seconds and Bytes 1.35G).  This could indicate a bug, or perhaps it's just a very complicated solution.

with(VectorCalculus); SetCoordinates(cartesian[x,y,z]);
F:= VectorField(<f1(x,y,z),f2(x,y,z),f3(x,y,z)>);
eq:= Laplacian(F) + 2*b*c*Curl(F) + c*F;
pdsolve(convert(eq,set));

However, you might look for solutions of the following form:

F = <a1,a2,a3>*exp(k*x)

(and those obtained by symmetry from that, and linear combinations...)

Maple does actually have a goto.  For example:

> F:= proc(x)
     if x = 2 then goto(b) end if;
     y:= 1;
     b:
     y
    end proc;

I think most of the goto's in Knuth can be replaced by judicious use of do loops together with break, next and return.  There are tricky cases, however, such as if you're in a triply nested loop and you want to break out of the two inner loops.

There isn't a case or switch statement in Maple.   Instead, you can use

if ... then ...
elif ... then ...

...

else ...

end if

1) Why can't you use the standard Bessel function?

2) You might try using RootFinding[NextZero].  For example:

> RootFinding[NextZero](x -> BesselJ(0,x), 0);

        2.404825557

One reason why garbage collection alone might not release enough memory is that large amounts of memory may be tied up in remember tables of procedures.  You might try sniffmem from my Maple Advisor Database,
http://www.math.ubc.ca/~israel/advisor.

The problem is that the left side of an assignment statement is not evaluated, so you can't really use seq on the left side.  However, you can use assign:

> seq(assign('A[i,i]'=i+2), i=1..3);

Let u = sqrt(1-2*x), v = sqrt(1+2*x), w = sqrt(1-4*x^2) = sqrt(u^2 * v^2) = (+/-) u*v.
Note that Maple uses the principal branch of the square root.  The branch cut for u is the real interval (1/2, +infinity), the branch cut for v is the real interval (-infinity, -1/2), and the branch cut for w is the union of those.  The complement of the union of the branch cuts is connected; since u*v = w rather than -w at some points of this set (e.g. x = 0), it will be w throughout this set.   So the only possibility to not have u*v = w is on the branch cuts.  For x in the interval (1/2, infinity), v is real and positive, and both u and w are on the positive imaginary axis, so u*v = w, not -w.  For x in the interval (-infinity, -1/2), u is real and positive, and both v and w are on the positive imaginary axis, so again u*v = w.  We conclude that under Maple's conventions, u*v = w on the whole complex plane.  Thus a/b = (2*x-1)/w * v/u = (2*x-1)/u^2 = -1 everywhere (excepting, of course, the points where it involves division by 0).