  Numerical bifurcation analysis, and in particular two-parameter continuation,
is used in consort with numerical simulation to reveal complicated dynamics in
the Mackey-Glass equation for moderate values of the delay close to the onset
of chaos. In particular a cusp bifurcation of periodic orbits and resulting
branches of folds of periodic orbits effectively partition the parameter space
into regions where different behaviours are seen. The cusp bifurcation leads
directly to bistability between periodic orbits, and subsequently to
bistability between a periodic orbit and a chaotic attractor. This leads to two
different mechanisms by which the chaotic attractor is destroyed in a global
bifurcation with a periodic orbit in either an interior crisis or a boundary
crisis. In another part of parameter space a sequence of subcritical
period-doublings is found to give rise to bistability between a periodic orbit
and a chaotic attractor. Torus bifurcations, and a codimension-two fold-flip
bifurcation are also identified, and Lyapunov exponent computations are used to
determine chaotic regions and attractor dimension.
