def METHOD_NAME(self, n=3, p=0, plot=True):
    pde = CosCosData()
    mesh = pde.init_mesh(n=n, meshtype='tri')
    space = RaviartThomasFiniteElementSpace2d(mesh, p=p)
    udof = space.number_of_global_dofs()
    pdof = space.smspace.number_of_global_dofs()
    gdof = udof + pdof
    uh = space.function()
    ph = space.smspace.function()
    A = space.stiff_matrix()
    B = space.div_matrix()
    F1 = space.smspace.source_vector(pde.source)
    AA = bmat([[A, -B], [-B.T, None]], format='csr')
    if True:
        F0 = -space.set_neumann_bc(pde.dirichlet)
        FF = np.r_['0', F0, F1]
        x = spsolve(AA, FF).reshape(-1)
        uh[:] = x[:udof]
        ph[:] = x[udof:]
        error0 = space.integralalg.L2_error(pde.flux, uh)
        def f(bc):
            xx = mesh.bc_to_point(bc)
            return (pde.solution(xx) - ph(xx))**2
        error1 = space.integralalg.integral(f)
        print(error0, error1)
    else:
        isBdDof = -space.set_dirichlet_bc(uh, pde.neumann)
        x = np.r_['0', uh, ph] 
        isBdDof = np.r_['0', isBdDof, np.zeros(pdof, dtype=np.bool_)]
        FF = np.r_['0', np.zeros(udof, dtype=np.float64), F1]
        FF -= AA@x
        bdIdx = np.zeros(gdof, dtype=np.int)
        bdIdx[isBdDof] = 1
        Tbd = spdiags(bdIdx, 0, gdof, gdof)
        T = spdiags(1-bdIdx, 0, gdof, gdof)
        AA = T@AA@T + Tbd
        FF[isBdDof] = x[isBdDof]
        x[:] = spsolve(AA, FF)
        uh[:] = x[:udof]
        ph[:] = x[udof:]
        error0 = space.integralalg.L2_error(pde.flux, uh)
        def f(bc):
            xx = mesh.bc_to_point(bc)
            return (pde.solution(xx) - ph(xx))**2
        error1 = space.integralalg.integral(f)
        print(error0, error1)
    if plot:
        box = [-0.5, 1.5, -0.5, 1.5]
        fig = plt.figure()
        axes = fig.gca()
        mesh.add_plot(axes, box=box)
        #mesh.find_node(axes, showindex=True)
        #mesh.find_edge(axes, showindex=True)
        #mesh.find_cell(axes, showindex=True)
        node = ps.reshape(-1, 2)
        uv = uh.reshape(-1, 2)
        axes.quiver(node[:, 0], node[:, 1], uv[:, 0], uv[:, 1])
        plt.show()