diff --git a/demos/cahnhilliard/demo_cahnhilliard.py.rst b/demos/cahnhilliard/demo_cahnhilliard.py.rst
index dc8d6c0..658efcf 100644
--- a/demos/cahnhilliard/demo_cahnhilliard.py.rst
+++ b/demos/cahnhilliard/demo_cahnhilliard.py.rst
@@ -34,6 +34,7 @@ For simplicity, we'll use a direct solver at each time step.
 Boilerplate imports::
 
   from firedrake import *
+  from firedrake.pyplot import tripcolor
   import numpy as np
   import os
   from irksome import Dt, GaussLegendre, MeshConstant, TimeStepper
diff --git a/demos/demo_nitsche_heat.py b/demos/demo_nitsche_heat.py
index eb37da2..dc3ad4b 100644
--- a/demos/demo_nitsche_heat.py
+++ b/demos/demo_nitsche_heat.py
@@ -23,7 +23,8 @@
 # MMS works on symbolic differentiation of true solution, not weak form
 rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
 
-u = interpolate(uexact, V)
+u = Function(V)
+u.interpolate(uexact)
 
 # define the variational form once outside the loop
 h = CellSize(msh)
diff --git a/demos/heat/demo_heat.py.rst b/demos/heat/demo_heat.py.rst
index 2df1e3e..60fe4d0 100644
--- a/demos/heat/demo_heat.py.rst
+++ b/demos/heat/demo_heat.py.rst
@@ -77,7 +77,8 @@ This defines the right-hand side using the method of manufactured solutions::
 We define the initial condition for the fully discrete problem, which
 will get overwritten at each time step::
 
-  u = interpolate(uexact, V)
+  u = Function(V)
+  u.interpolate(uexact)
 
 Now, we will define the semidiscrete variational problem using
 standard UFL notation, augmented by the ``Dt`` operator from Irksome::
diff --git a/demos/heat/demo_heat_dirk.py.rst b/demos/heat/demo_heat_dirk.py.rst
index fd1050c..78273ce 100644
--- a/demos/heat/demo_heat_dirk.py.rst
+++ b/demos/heat/demo_heat_dirk.py.rst
@@ -76,7 +76,8 @@ This defines the right-hand side using the method of manufactured solutions::
 We define the initial condition for the fully discrete problem, which
 will get overwritten at each time step::
 
-  u = interpolate(uexact, V)
+  u = Function(V)
+  u.interpolate(uexact)
 
 Now, we will define the semidiscrete variational problem using
 standard UFL notation, augmented by the ``Dt`` operator from Irksome::
diff --git a/demos/lowlevel/demo_lowlevel_homogbc.py.rst b/demos/lowlevel/demo_lowlevel_homogbc.py.rst
index e80ba12..11b4e70 100644
--- a/demos/lowlevel/demo_lowlevel_homogbc.py.rst
+++ b/demos/lowlevel/demo_lowlevel_homogbc.py.rst
@@ -46,7 +46,8 @@ Continuing::
 
   rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
 
-  u = interpolate(uexact, V)
+  u = Function(V)
+  u.interpolate(uexact)
 
   v = TestFunction(V)
   F = inner(Dt(u), v)*dx + inner(grad(u), grad(v))*dx - inner(rhs, v)*dx
diff --git a/demos/lowlevel/demo_lowlevel_inhomogbc.py.rst b/demos/lowlevel/demo_lowlevel_inhomogbc.py.rst
index b7791e6..509e235 100644
--- a/demos/lowlevel/demo_lowlevel_inhomogbc.py.rst
+++ b/demos/lowlevel/demo_lowlevel_inhomogbc.py.rst
@@ -28,7 +28,8 @@ Imports::
   uexact = exp(-t) * cos(pi * x) * sin(pi * y)
   rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
 
-  u = interpolate(uexact, V)
+  u = Function(V)
+  u.interpolate(uexact)
 
   v = TestFunction(V)
   F = inner(Dt(u), v)*dx + inner(grad(u), grad(v))*dx - inner(rhs, v)*dx
diff --git a/demos/navier_stokes/navier_stokes_demo.py b/demos/navier_stokes/navier_stokes_demo.py
deleted file mode 100644
index 314d722..0000000
--- a/demos/navier_stokes/navier_stokes_demo.py
+++ /dev/null
@@ -1,129 +0,0 @@
-from firedrake import *  # noqa: F403
-from irksome import LobattoIIIC, getForm, Dt, MeshConstant, TimeStepper
-import numpy
-
-ButcherTableau = LobattoIIIC(3)
-
-# Corners of the box.
-x0 = 0.0
-x1 = 2.2
-y0 = 0.0
-y1 = 0.41
-
-#Mesh
-stepfile = "circle.step"
-# Close to what Turek's time-dependent benchmark has on coarsest mesh
-h = 0.05 
-order = 2
-lvl = 1
-mh = OpenCascadeMeshHierarchy(stepfile, element_size=h,
-                              levels=lvl, order=order, cache=False,
-                              verbose=True)
-
-#Space
-msh = mh[-1]
-
-
-V = VectorFunctionSpace(msh, "CG", 4)
-W = FunctionSpace(msh, "DG", 3)
-Z = V * W
-
-v, w = TestFunctions(Z)
-x, y = SpatialCoordinate(msh)
-
-up = Function(Z)
-u, p = split(up)
-
-n= FacetNormal(msh)
-
-MC = MeshConstant(msh)
-dt = MC.Constant(1.0/1600)
-t = MC.Constant(0.0)
-nu = Constant(0.001)
-rho=1.0
-r=0.05
-Umean=1.0
-L=0.1
-
-# define the variational form once outside the loop
-gamma = Constant(1.e2)
-F = (inner(Dt(u), v) * dx
-     +inner(dot(grad(u), u), v) * dx
-     #+ gamma * inner(cell_avg(div(u)), cell_avg(div(v))) * dx
-     + nu * inner(grad(u), grad(v)) * dx
-     - inner(p, div(v)) * dx
-     + inner(div(u), w) * dx)
-
-# boundary conditions are specified for each subspace
-UU = 1.5 * sin(pi * t / 8)
-bcs = [DirichletBC(Z.sub(0),
-                   as_vector([4 * UU * y * (y1 - y) / (y1**2), 0]), (4,)),
-       DirichletBC(Z.sub(0), Constant((0, 0)), (1, 3, 5))]
-
-#nullspace = MixedVectorSpaceBasis(
-#    Z, [Z.sub(0), VectorSpaceBasis(constant=True)])
-
-s_param={'ksp_type': 'preonly',
-         'pc_type': 'lu',
-         'pc_factor_shift_type': 'inblocks',
-         'mat_type': 'aij',
-         #'snes_monitor': None
-         }
-
-
-
-# the TimeStepper object builds the UFL for the multi-stage RK method
-# and sets up a solver.  It also provides a method for updating the solution
-# at each time step.
-# It overwites the UFL constant t with the current time value at each
-# step.  u is also updated in place.
-stepper = TimeStepper(F, ButcherTableau, t, dt, up, bcs=bcs,
-                      solver_parameters=s_param)
-
-times = []
-CD = []
-CL = []
-
-ut1 = dot(u, as_vector([n[1], -n[0]]))
-
-while (float(t) < 8):
-    # Update step
-    stepper.advance()
-    print(float(t))
-    t.assign(float(t) + float(dt))
-
-    # Numerical quantities
-    FD= assemble((nu*dot(grad(ut1), n)*n[1] - p*n[0])*ds(5) )
-    FL= assemble(-(nu*inner(grad(ut1), n)*n[0]+p*n[1])*ds(5) )
-
-    CD+= [-2 * FD / (Umean**2) / L] #Drag coefficient
-    CL+= [-2 * FL / (Umean**2) / L] #Lift coefficient
-    times.append(float(t))
-
-
-times = numpy.asarray(times)
-CD = numpy.asarray(CD)
-CL = numpy.asarray(CL)
-
-iCD = numpy.argmax(CD)
-tCDmax = times[iCD]
-CDmax = CD[iCD]
-
-iCL = numpy.argmax(CL)
-tCLmax = times[iCL]
-CLmax = CL[iCL]
-
-#p_func=interpolate(p, W)
-#p1=p_func( [0.15,0.2] )
-#p2=p_func( [0.25,0.2] )
-#pdiff= p1-p2
-
-#Results
-print("Level", lvl)
-print("CD_max", CDmax)
-print("t(CD_max)", tCDmax)
-print("CL_max", CLmax)
-print("t(CL_max)", tCLmax)
-
-#u, p = up.split()
-#print(p((0.15, .2)) - p((0.250, 0.2)))
diff --git a/demos/navier_stokes/nse_steady_demo.py b/demos/navier_stokes/nse_steady_demo.py
deleted file mode 100644
index 43120ac..0000000
--- a/demos/navier_stokes/nse_steady_demo.py
+++ /dev/null
@@ -1,130 +0,0 @@
-from firedrake import *  # noqa: F403
-import numpy
-
-
-class DGMassInv(PCBase):
-
-    def initialize(self, pc):
-        _, P = pc.getOperators()
-        appctx = self.get_appctx(pc)
-        V = dmhooks.get_function_space(pc.getDM())
-        # get function spaces
-        u = TrialFunction(V)
-        v = TestFunction(V)
-        massinv = assemble(Tensor(inner(u, v)*dx).inv)
-        self.massinv = massinv.petscmat
-        self.nu = appctx["nu"]
-        self.gamma = appctx["gamma"]
-
-    def update(self, pc):
-        pass
-
-    def apply(self, pc, x, y):
-        self.massinv.mult(x, y)
-        scaling = float(self.nu) + float(self.gamma)
-        y.scale(-scaling)
-
-    def applyTranspose(self, pc, x, y):
-        raise NotImplementedError("Sorry!")
-
-
-# Corners of the box.
-x0 = 0.0
-x1 = 2.2
-y0 = 0.0
-y1 = 0.41
-
-#Mesh
-stepfile = "circle.step"
-# Close to what Turek's time-dependent benchmark has on coarsest mesh
-h = 0.1 
-order = 2
-lvl = 2
-mh = OpenCascadeMeshHierarchy(stepfile, element_size=h,
-                              levels=lvl, order=order, cache=False,
-                              verbose=True)
-
-msh = mh[-1]
-
-
-
-V = VectorFunctionSpace(msh, "CG", 4)
-W = FunctionSpace(msh, "DG", 3)
-Z = V * W
-
-v, w = TestFunctions(Z)
-x, y = SpatialCoordinate(msh)
-
-up = Function(Z)
-u, p = split(up)
-
-n= FacetNormal(msh)
-
-nu = Constant(0.001)
-rho=1.0
-r=0.05
-Umean=0.2
-L=0.1
-
-# define the variational form once outside the loop
-gamma = Constant(1.e2)
-F = (inner(dot(grad(u), u), v) * dx
-     + gamma * inner(div(u),div(v)) * dx
-     + nu * inner(grad(u), grad(v)) * dx
-     - inner(p, div(v)) * dx
-     + inner(div(u), w) * dx)
-
-# boundary conditions are specified for each subspace
-UU = Constant(0.3)
-bcs = [DirichletBC(Z.sub(0),
-                   as_vector([4 * UU * y * (y1 - y) / (y1**2), 0]), (4,)),
-       DirichletBC(Z.sub(0), Constant((0, 0)), (1, 3, 5))]
-
-
-s_param = {'snes_monitor': None,
-           'ksp_type': 'fgmres',
-           'ksp_monitor': None,
-           'pc_type': 'fieldsplit',
-           'pc_fieldsplit_type': 'schur',
-           'pc_fieldsplit_schur_fact_type': 'full',
-           'fieldsplit_0': {
-               'ksp_type': 'preonly',
-               'pc_type': 'lu'},
-           'fieldsplit_1': {
-               'ksp_type': 'gmres',
-               'ksp_max_it': 1,
-               'pc_type': 'python',
-               'pc_python_type': '__main__.DGMassInv'
-               }}
-
-
-appctx = {'nu': nu, 'gamma': float(gamma)}
-ut1 = dot(u, as_vector([n[1], -n[0]]))
-
-solve(F==0, up, bcs=bcs, solver_parameters=s_param, appctx=appctx)
-
-FD= assemble((nu*dot(grad(ut1), n)*n[1] - p*n[0])*ds(5) )
-FL= assemble(-(nu*inner(grad(ut1), n)*n[0]+p*n[1])*ds(5) )
-
-CD = [-2 * FD / (Umean**2) / L] #Drag coefficient
-CL = [-2 * FL / (Umean**2) / L] #Lift coefficient
-
-print("Level", lvl)
-print("CD", CD)
-print("CL", CL)
-
-u, p = up.split()
-
-import time
-tstart = time.time()
-pdiffat = p.at((0.15, .2)) - p.at((0.250, 0.2))
-tat = time.time() - tstart
-print("time using at ", tat)
-tstart = time.time()
-pdiffnotat = p((0.15, .2)) - p((0.25, .2))
-tnotat = time.time() - tstart
-print("time without at ", tnotat)
-
-
-# u, p = up.split()
-# File("nse_steady.pvd").write(u, p)
diff --git a/demos/preconditioning/demo_heat_mg.py.rst b/demos/preconditioning/demo_heat_mg.py.rst
index 8d0b61e..24d844d 100644
--- a/demos/preconditioning/demo_heat_mg.py.rst
+++ b/demos/preconditioning/demo_heat_mg.py.rst
@@ -57,7 +57,8 @@ are just as for the regular heat equation demo::
   uexact = B * atan(t)*(pi / 2.0 - atan(S * (R - t)))
   rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
 
-  u = interpolate(uexact, V)
+  u = Function(V)
+  u.interpolate(uexact)
   v = TestFunction(V)
   F = inner(Dt(u), v)*dx + inner(grad(u), grad(v))*dx - inner(rhs, v)*dx
   bc = DirichletBC(V, 0, "on_boundary")
diff --git a/demos/preconditioning/demo_heat_pc.py.rst b/demos/preconditioning/demo_heat_pc.py.rst
index 82711ff..9145011 100644
--- a/demos/preconditioning/demo_heat_pc.py.rst
+++ b/demos/preconditioning/demo_heat_pc.py.rst
@@ -69,8 +69,8 @@ Common set-up for the problem::
 
   uexact = B * atan(t)*(pi / 2.0 - atan(S * (R - t)))
   rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
-  u = interpolate(uexact, V)
-
+  u = Function(V)
+  u.interpolate(uexact)
   v = TestFunction(V)
   F = inner(Dt(u), v)*dx + inner(grad(u), grad(v))*dx - inner(rhs, v)*dx
 
diff --git a/demos/preconditioning/demo_heat_rana.py.rst b/demos/preconditioning/demo_heat_rana.py.rst
index b95a04a..c5b0a21 100644
--- a/demos/preconditioning/demo_heat_rana.py.rst
+++ b/demos/preconditioning/demo_heat_rana.py.rst
@@ -68,7 +68,8 @@ Common set-up for the problem::
 
   uexact = B * atan(t)*(pi / 2.0 - atan(S * (R - t)))
   rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
-  u = interpolate(uexact, V)
+  u = Function(V)
+  u.interpolate(uexact)
 
   v = TestFunction(V)
   F = inner(Dt(u), v)*dx + inner(grad(u), grad(v))*dx - inner(rhs, v) * dx
diff --git a/irksome/dirk_stepper.py b/irksome/dirk_stepper.py
index ddd1822..0527612 100644
--- a/irksome/dirk_stepper.py
+++ b/irksome/dirk_stepper.py
@@ -2,7 +2,8 @@
 from firedrake import Constant, DirichletBC, Function
 from firedrake import NonlinearVariationalProblem as NLVP
 from firedrake import NonlinearVariationalSolver as NLVS
-from firedrake import interpolate, split
+from firedrake import assemble, split, project
+from firedrake.__future__ import interpolate
 from ufl.constantvalue import as_ufl
 
 from .deriv import TimeDerivative
@@ -110,10 +111,10 @@ def getFormDIRK(F, butch, t, dt, u0, bcs=None):
         bcarg = as_ufl(bc._original_arg)
         bcarg_stage = replace(bcarg, {t: t+c*dt})
         try:
-            gdat = interpolate(bcarg, Vbc)
+            gdat = assemble(interpolate(bcarg, Vbc))
             gmethod = lambda gd, gc: gd.interpolate(gc)
         except:  # noqa: E722
-            gdat = interpolate(bcarg, Vbc)
+            gdat = assemble(project(bcarg, Vbc))
             gmethod = lambda gd, gc: gd.project(gc)
 
         new_bc = DirichletBC(Vbc, gdat, bc.sub_domain)
diff --git a/irksome/getForm.py b/irksome/getForm.py
index 72c8e92..69d6625 100644
--- a/irksome/getForm.py
+++ b/irksome/getForm.py
@@ -3,7 +3,8 @@
 
 import numpy
 from firedrake import (Constant, DirichletBC, Function, TestFunction,
-                       interpolate, project, split)
+                       assemble, project, split)
+from firedrake.__future__ import interpolate
 from ufl import diff
 from ufl.algorithms import expand_derivatives
 from ufl.classes import Zero
@@ -18,7 +19,7 @@ def __init__(self, V, gcur, u0, u0_mult, i, t, dt):
             if V.parent.index is None:  # but not part of a MFS
                 sub = V.component
                 try:
-                    gdat = interpolate(gcur-u0_mult[i]*u0.sub(sub), V)
+                    gdat = assemble(interpolate(gcur-u0_mult[i]*u0.sub(sub), V))
                     gmethod = lambda g, u: gdat.interpolate(g-u0_mult[i]*u.sub(sub))
                 except:  # noqa: E722
                     gdat = project(gcur-u0_mult[i]*u0.sub(sub), V)
@@ -27,7 +28,7 @@ def __init__(self, V, gcur, u0, u0_mult, i, t, dt):
                 sub0 = V.parent.index
                 sub1 = V.component
                 try:
-                    gdat = interpolate(gcur-u0_mult[i]*u0.sub(sub0).sub(sub1), V)
+                    gdat = assemble(interpolate(gcur-u0_mult[i]*u0.sub(sub0).sub(sub1), V))
                     gmethod = lambda g, u: gdat.interpolate(g-u0_mult[i]*u.sub(sub0).sub(sub1))
                 except:  # noqa: E722
                     gdat = project(gcur-u0_mult[i]*u0.sub(sub0).sub(sub1), V)
@@ -35,7 +36,7 @@ def __init__(self, V, gcur, u0, u0_mult, i, t, dt):
         else:  # V is not a bit of a VFS
             if V.index is None:  # not part of MFS, either
                 try:
-                    gdat = interpolate(gcur-u0_mult[i]*u0, V)
+                    gdat = assemble(interpolate(gcur-u0_mult[i]*u0, V))
                     gmethod = lambda g, u: gdat.interpolate(g-u0_mult[i]*u)
                 except:  # noqa: E722
                     gdat = project(gcur-u0_mult[i]*u0, V)
@@ -43,7 +44,7 @@ def __init__(self, V, gcur, u0, u0_mult, i, t, dt):
             else:  # part of MFS
                 sub = V.index
                 try:
-                    gdat = interpolate(gcur-u0_mult[i]*u0.sub(sub), V)
+                    gdat = assemble(interpolate(gcur-u0_mult[i]*u0.sub(sub), V))
                     gmethod = lambda g, u: gdat.interpolate(g-u0_mult[i]*u.sub(sub))
                 except:  # noqa: E722
                     gdat = project(gcur-u0_mult[i]*u0.sub(sub), V)
diff --git a/irksome/stage.py b/irksome/stage.py
index 4c33881..35b869f 100644
--- a/irksome/stage.py
+++ b/irksome/stage.py
@@ -5,7 +5,8 @@
 import numpy as np
 from firedrake import (Constant, DirichletBC, Function,
                        NonlinearVariationalProblem, NonlinearVariationalSolver,
-                       TestFunction, dx, inner, interpolate, project, split)
+                       TestFunction, assemble, dx, inner, project, split)
+from firedrake.__future__ import interpolate
 from numpy import vectorize
 from ufl.classes import Zero
 from ufl.constantvalue import as_ufl
@@ -242,7 +243,7 @@ def getFormStage(F, butch, u0, t, dt, bcs=None, splitting=None,
         bcarg = as_ufl(bc._original_arg)
         for i in range(num_stages):
             try:
-                gdat = interpolate(bcarg, Vsp)
+                gdat = assemble(interpolate(bcarg, Vsp))
                 gmethod = lambda gd, gc: gd.interpolate(gc)
             except:  # noqa: E722
                 gdat = project(bcarg, Vsp)
@@ -292,7 +293,7 @@ def getFormStage(F, butch, u0, t, dt, bcs=None, splitting=None,
 
         bcarg = as_ufl(bc._original_arg)
         try:
-            gdat = interpolate(bcarg, Vsp)
+            gdat = assemble(interpolate(bcarg, Vsp))
             gmethod = lambda gd, gc: gd.interpolate(gc)
         except:  # noqa: E722
             gdat = project(bcarg, Vsp)
diff --git a/tests/test_curl.py b/tests/test_curl.py
index 2e2c09c..bb79310 100644
--- a/tests/test_curl.py
+++ b/tests/test_curl.py
@@ -1,5 +1,6 @@
 import pytest
 from firedrake import *
+from firedrake.__future__ import interpolate
 from irksome import GaussLegendre, Dt, MeshConstant, TimeStepper
 from irksome.tools import AI, IA
 from ufl.algorithms.ad import expand_derivatives
@@ -26,7 +27,7 @@ def curltest(N, deg, butcher_tableau, splitting):
     uexact = as_vector([t + 2*t*x + 4*t*y + 3*t*(y**2) + 2*t*x*y, 7*t + 5*t*x + 6*t*y - 3*t*x*y - 2*t*(x**2)])
     rhs = expand_derivatives(diff(uexact, t)) + curl(curl(uexact))
 
-    u = interpolate(uexact, V)
+    u = assemble(interpolate(uexact, V))
 
     v = TestFunction(V)
     n = FacetNormal(msh)
diff --git a/tests/test_dirichletbc.py b/tests/test_dirichletbc.py
index dbc0fa6..29e4935 100644
--- a/tests/test_dirichletbc.py
+++ b/tests/test_dirichletbc.py
@@ -35,7 +35,9 @@ def test_1d_heat_dirichletbc(butcher_tableau, stage_type):
         + ((x - x0) / x1) * (u_1 - u_0)
     )
     rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
-    u = interpolate(uexact, V)
+    u = Function(V)
+    u.interpolate(uexact)
+
     v = TestFunction(V)
     F = (
         inner(Dt(u), v) * dx
diff --git a/tests/test_dirk.py b/tests/test_dirk.py
index b21c865..5be278d 100644
--- a/tests/test_dirk.py
+++ b/tests/test_dirk.py
@@ -34,7 +34,8 @@ def test_1d_heat_dirichletbc(butcher_tableau):
         + ((x - x0) / x1) * (u_1 - u_0)
     )
     rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
-    u = interpolate(uexact, V)
+    u = Function(V)
+    u.interpolate(uexact)
     v = TestFunction(V)
     F = (
         inner(Dt(u), v) * dx
@@ -78,8 +79,11 @@ def test_1d_heat_neumannbc(butcher_tableau):
 
     uexact = cos(pi*x)*exp(-(pi**2)*t)
     rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
-    u_dirk = interpolate(uexact, V)
-    u = interpolate(uexact, V)
+    u_dirk = Function(V)
+    u = Function(V)
+    u_dirk.interpolate(uexact)
+    u.interpolate(uexact)
+
     v = TestFunction(V)
     F = (
         inner(Dt(u), v) * dx
@@ -120,8 +124,10 @@ def test_1d_heat_homogdbc(butcher_tableau):
 
     uexact = sin(pi*x)*exp(-(pi**2)*t)
     rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
-    u_dirk = interpolate(uexact, V)
-    u = interpolate(uexact, V)
+    u_dirk = Function(V)
+    u = Function(V)
+    u_dirk.interpolate(uexact)
+    u.interpolate(uexact)
     v = TestFunction(V)
     F = (
         inner(Dt(u), v) * dx
@@ -167,8 +173,10 @@ def test_1d_vectorheat_componentBC(butcher_tableau):
     uexact = as_vector([sin(pi*x/2)*exp(-(pi**2)*t/4),
                         cos(pi*x/2)*exp(-(pi**2)*t/4)])
     rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
-    u_dirk = interpolate(uexact, V)
-    u = interpolate(uexact, V)
+    u_dirk = Function(V)
+    u = Function(V)
+    u_dirk.interpolate(uexact)
+    u.interpolate(uexact)
     v = TestFunction(V)
     F = (
         inner(Dt(u), v) * dx
diff --git a/tests/test_inhomogbc.py b/tests/test_inhomogbc.py
index 5a5adce..e1e2443 100644
--- a/tests/test_inhomogbc.py
+++ b/tests/test_inhomogbc.py
@@ -18,7 +18,8 @@ def heat_inhomog(N, deg, butcher_tableau, stage_type, splitting):
     uexact = t*(x+y)
     rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
 
-    u = interpolate(uexact, V)
+    u = Function(V)
+    u.interpolate(uexact)
 
     v = TestFunction(V)
     F = inner(Dt(u), v)*dx + inner(grad(u), grad(v))*dx - inner(rhs, v)*dx
diff --git a/tests/test_pc.py b/tests/test_pc.py
index cc9a564..b4ef297 100644
--- a/tests/test_pc.py
+++ b/tests/test_pc.py
@@ -1,8 +1,8 @@
 import numpy
 import pytest
-from firedrake import (DirichletBC, FunctionSpace, SpatialCoordinate,
+from firedrake import (DirichletBC, Function, FunctionSpace, SpatialCoordinate,
                        TestFunction, UnitSquareMesh, diff, div, dx, errornorm,
-                       grad, inner, interpolate)
+                       grad, inner)
 from irksome import Dt, MeshConstant, LobattoIIIC, RadauIIA, TimeStepper
 from ufl.algorithms.ad import expand_derivatives
 
@@ -10,7 +10,8 @@
 
 
 def Fubc(V, uexact, rhs):
-    u = interpolate(uexact, V)
+    u = Function(V)
+    u.interpolate(uexact)
     v = TestFunction(V)
     F = inner(Dt(u), v)*dx + inner(grad(u), grad(v))*dx - inner(rhs, v)*dx
 
diff --git a/tests/test_split.py b/tests/test_split.py
index d1ae834..679a7d7 100644
--- a/tests/test_split.py
+++ b/tests/test_split.py
@@ -4,7 +4,7 @@
                        TestFunction, TrialFunctions, UnitIntervalMesh,
                        UnitSquareMesh, VectorElement, VectorSpaceBasis, action,
                        as_vector, assemble, derivative, div, dot, dx,
-                       errornorm, grad, inner, interpolate, pi, sin, split)
+                       errornorm, grad, inner, pi, sin, split)
 from irksome import Dt, MeshConstant, RadauIIA, TimeStepper
 from irksome.tools import AI, IA
 
@@ -26,8 +26,10 @@ def test_diffreact(splitting):
 
     uIC = 2 + sin(2 * pi * x)
 
-    u_imp = interpolate(uIC, V)
-    u_split = interpolate(uIC, V)
+    u_imp = Function(V)
+    u_split = Function(V)
+    u_imp.interpolate(uIC)
+    u_split.interpolate(uIC)
 
     bcs = DirichletBC(V, Constant(2), "on_boundary")
 
diff --git a/tests/test_subdomainbc.py b/tests/test_subdomainbc.py
index be4486a..cca83cc 100644
--- a/tests/test_subdomainbc.py
+++ b/tests/test_subdomainbc.py
@@ -18,7 +18,8 @@ def heat_subdomainbc(N, deg, butcher_tableau, splitting=AI):
     uexact = t*(x+y)
     rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
 
-    u = interpolate(uexact, V)
+    u = Function(V)
+    u.interpolate(uexact)
 
     v = TestFunction(V)
     n = FacetNormal(msh)
@@ -56,7 +57,8 @@ def heat_componentbc(N, deg, butcher_tableau, splitting=AI):
     uexact = as_vector([(x**2-2*x)*t, (1-x**2)*t])
     rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
 
-    u = interpolate(uexact, V)
+    u = Function(V)
+    u.interpolate(uexact)
 
     v = TestFunction(V)
     n = FacetNormal(msh)