def METHOD_NAME(p, n, maxit):
    from fealpy.pde.elliptic_1d import SinPDEData as PDE
    from fealpy.mesh import IntervalMesh
    pde = PDE()
    domain = pde.domain()
    em = np.zeros((2, maxit), dtype=np.float64)
    for i in range(maxit):
        mesh = IntervalMesh.from_interval_domain(domain, nx=n * 2**i)
        space = Space(mesh, p=p)
        bform = BilinearForm(space)
        bform.add_domain_integrator(ScalarDiffusionIntegrator())
        bform.assembly()
        lform = LinearForm(space)
        lform.add_domain_integrator(ScalarSourceIntegrator(pde.source, q=p+3))
        lform.assembly()
        A = bform.get_matrix()
        f = lform.get_vector() 
        ospace = OldSpace(mesh, p=p)
        OA = ospace.stiff_matrix()
        Of = ospace.source_vector(pde.source, q=p+3)
        np.testing.assert_allclose(A.toarray(), OA.toarray())
        #np.testing.assert_allclose(f, Of)
        bc = DirichletBC(space, pde.dirichlet)
        uh = space.function()
        A, f = bc.apply(A, f, uh)
        uh[:] = spsolve(A, f)
        em[0, i] = mesh.error(pde.solution, uh, q=p+3)
        em[1, i] = mesh.error(pde.gradient, uh.grad_value, q=p+3)
    ratio = em[:, 0:-1]/em[:, 1:]
    print(em)
    print(ratio)
    assert np.abs(ratio[0, -1] - 2**(p+1)) < 0.1
    assert np.abs(ratio[1, -1] - 2**p) < 0.1