Matplotlib fill area between contour lines where one contour line is made of two disjoint curves

Question

I have a function f and I would like to color-fill between two contours. One contour is shown in blue in the figure on the left. The contour consists in two disjoint curves. Another contour is shown in green in the middle figure, and consists in a single curve. Both contours are shown together in the final image on the right.

How can I fill in between these two contours (generated using ax.contour()) when one contour is made of two disjoin curves?

Minimal Working Example

The following code creates the function of interest, generates a grid of values, and then plots and stores the contours.

import numpy as np
from numpy import apply_along_axis, meshgrid, arange, vstack
import matplotlib.pyplot as plt

# Function whose level-sets we are interested in
f = lambda x: x[1]**2 + 3*(x[0]**2)*(x[0]**2 - 1)
f_broadcasted = lambda x: apply_along_axis(f, 1, x)

# Generate grid of (x, y, z) values to plot contour
xlims, ylims = [-2, 2], [-2, 2]
x, y = meshgrid(arange(*xlims, 0.005), arange(*ylims, 0.005))
xys = vstack((x.flatten(), y.flatten())).T
z = f_broadcasted(xys).reshape(x.shape)

# Plot contours
fig, ax = plt.subplots()
contours = ax.contour(x, y, z, levels=[-0.02, 0.02], colors=['forestgreen', 'royalblue'])
plt.show()

The output of this is basically similar to the right-most figure above. I tried using

ax.fill_between(*contours[1][0].T, contours[0][1])

but the problem is that they have different sizes.

Basically my desired output could be achieved (in an ugly way) by plotting MANY contours

fig, ax = plt.subplots()
contours = ax.contour(x, y, z, levels=np.linspace(-0.02, 0.02, num=1000), colors='black', linestyles='-')
plt.show()

but the problem is that this takes quite a while and the final product is not particularly neat.

Answer 1

Using contourf:

import numpy as np
from numpy import apply_along_axis, meshgrid, arange, vstack
import matplotlib.pyplot as plt

# Function whose level-sets we are interested in
f = lambda x: x[1]**2 + 3*(x[0]**2)*(x[0]**2 - 1)
f_broadcasted = lambda x: apply_along_axis(f, 1, x)

# Generate grid of (x, y, z) values to plot contour
xlims, ylims = [-2, 2], [-2, 2]
x, y = meshgrid(arange(*xlims, 0.005), arange(*ylims, 0.005))
xys = vstack((x.flatten(), y.flatten())).T
z = f_broadcasted(xys).reshape(x.shape)

# Plot contours
fig, ax = plt.subplots()
contours = ax.contourf(x, y, z, levels=[-0.02, 0.02], colors='forestgreen')
plt.show()

Output:

Combining contour and contourf for a nicer display:

contours = ax.contourf(x, y, z, levels=[-0.02, 0.02], colors='forestgreen')
ax.contour(x, y, z, levels=[-0.02, 0.02],
           colors=['forestgreen', 'royalblue'], linestyles=['-'])

Output:

Using the nchunk parameter for better antialiasing, still with some artifacts:

contours = ax.contourf(x, y, z, levels=[-0.02, 0.02],
                       antialiased=True, nchunk=20,
                       colors='forestgreen')

Thanks to @JohanC for the suggestion.