plotting modes from demo_waveguide.py example

I have trouble reproducing the waveguide cutoff modes plots using the demo provided. A simple plot(f) command does not produce anything resembling a mode although the calculated cutoff frequency is the same as one in Fenics book. The code is given below.

    from dolfin import *

    # Test for PETSc and SLEPc
    if not has_linear_algebra_backend("PETSc"):
        print "DOLFIN has not been configured with PETSc. Exiting."
        exit()

    if not has_slepc():
        print "DOLFIN has not been configured with SLEPc. Exiting."
        exit()

    # Make sure we use the PETSc backend
    parameters["linear_algebra_backend"] = "PETSc"

    # Create mesh
    width = 1.0
    height = 0.5
    mesh = RectangleMesh(0, 0, width, height, 4, 2)

    # Define the function space
    V = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 3)

    # Define the test and trial functions
    v = TestFunction(V)
    u = TrialFunction(V)

    # Define the forms - generates an generalized eigenproblem of the form
    # [S]{h} = k_o^2[T]{h}
    def curl_t(w):
        return Dx(w[1], 0) - Dx(w[0], 1)
    s = curl_t(v)*curl_t(u)*dx
    t = inner(v, u)*dx

    # Assemble the stiffness matrix (S) and mass matrix (T)
    S = PETScMatrix()
    T = PETScMatrix()
    assemble(s, tensor=S)
    assemble(t, tensor=T)

    # Solve the eigensystem
    esolver = SLEPcEigenSolver(S, T)
    esolver.parameters["spectrum"] = "smallest real"
    esolver.parameters["solver"] = "lapack"
    esolver.solve()

    cutoff = None
    for i in range(S.size(1)):
        (lr, lc) = esolver.get_eigenvalue(i)
        if lr > 1 and lc == 0:
        #extracting the eigenvector
        lr, lc, evec_r, evec_c = esolver.get_eigenpair(i)
            cutoff = sqrt(lr)
            break

    if cutoff is None:
        print "Unable to find dominant mode"
    else:
        print "Cutoff frequency:", cutoff

#simple try at plotting?

f = Function(V,evec_r)
plot(f,mode="color")

Another method would be to extract the components manually, but here I also got stuck. Here is my attempt.

import numpy as np                                                             
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt 

n = V.dim()
d = mesh.geometry().dim()

dof_coordinates = V.dofmap().tabulate_all_coordinates(mesh)                     
dof_coordinates.resize((n, d))                                                  
dof_x = dof_coordinates[:, 0]                                                   
dof_y = dof_coordinates[:, 1] 

fig = plt.figure()                                                              
ax = fig.add_subplot(111, projection='3d')                                      
ax.scatter(dof_x, dof_y, f.vector().array(), c='b', marker='.')                 
plt.show()

For Nedelec elements what do the dof_coordinates correspond to? Clearly these are not vertices. How would I convert the entries of the eigenvector (evec_r,evec_c) to the corresponding Ex and Ey fields? And finally, are these located at the vertices of at the mesh or at points specified by dof_coordinates?
Both of these methods do not work for me, and do not understand why. Any thoughts?