diff --git a/doc/user_other_numerical_integration.rst b/doc/user_other_numerical_integration.rst
index 3db3630d..176c2509 100644
--- a/doc/user_other_numerical_integration.rst
+++ b/doc/user_other_numerical_integration.rst
@@ -212,11 +212,8 @@ The usage pattern is as follows:
     rsq = grid.points[:,0]**2 + grid.points[:,1]**2 + grid.points[:,2]**2
     pot = np.sqrt(1 + rsq)
 
-    # Allocate an output array for the operator
-    fock = lf.create_two_index()
-
     # Actual computation
-    obasis.compute_grid_density_fock(grid.points, grid.weights, pot, fock)
+    fock = obasis.compute_grid_density_fock(grid.points, grid.weights, pot)
 
 
 Other chain rules are also implemented:
diff --git a/setup.py b/setup.py
index 8f959579..9e68d0ac 100755
--- a/setup.py
+++ b/setup.py
@@ -442,7 +442,6 @@ def lib_config_magic(prefix, libname, static_config={}, known_include_dirs=[]):
         'horton.gbasis', 'horton.gbasis.test',
         'horton.grid', 'horton.grid.test',
         'horton.io', 'horton.io.test',
-        'horton.matrix', 'horton.matrix.test',
         'horton.meanfield', 'horton.meanfield.test',
         'horton.part', 'horton.part.test',
         'horton.scripts', 'horton.scripts.test',
@@ -467,7 +466,6 @@ def lib_config_magic(prefix, libname, static_config={}, known_include_dirs=[]):
         ('include/horton/grid', glob('horton/grid/*.h')),
         ('include/horton/gbasis', glob('horton/gbasis/*.h')),
         ('include/horton/espfit', glob('horton/espfit/*.h')),
-        ('include/horton/matrix', glob('horton/matrix/*.h')),
     ],
     package_data={
         'horton': ['*.pxd'],