diff --git a/README.md b/README.md
index c142b8e..c8d3088 100644
--- a/README.md
+++ b/README.md
@@ -3,7 +3,7 @@
 A package for computing unassembled matrices appearing in $`C^0`$ continuous Galerkin
 finite element models.
 
-# Methods
+# Matrices in one coordinate (1D)
 
 This package computes matrices of the following form.
 ```math
@@ -45,6 +45,10 @@ N_{ijk} = \int^{x_u}_{x_l} P_i(x) Q_j(x) R_k(x) \rho(x) d x = \int^{1}_{-1} S_i(
 where $`R_i(x)`$ may be either of $`\Phi(x)`$ or $`d\Phi/dx`$ and
 $`U_i(z)`$ may be either of $`l_i(z)`$ or $`(1/s)d l_i(z)/d z`$.
 
+# p-adaptive Gaussian quadrature
+
+In order to guarantee that the matrices $`A_{ij}`$ with non-polynomial kernels $`\rho(x)`$ are calculated with sufficient accuracy, we use p-adaptive quadrature with Gauss-Legendre quadrature rules from `FastGaussQuadrature` at each iteration. We iteratively calculate the matrices with an increasingly large number of quadrature points and iterate until the result is constant up to absolute and relative tolerances. The inputs used to control this process are described in the section on usage below. For polynomial kernels we can deduce the number of quadrature points for exact results, avoiding the need for the iterative calculation in the pure polynomial case.
+
 # Usage
 
 To create these matrices, the user must supply the key grid
@@ -75,6 +79,16 @@ using FiniteElementMatrices
 coordinate = ElementCoordinates(x,scale,shift)
 M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
 ```
+Note that the order in which the `@enum` Lagrange polynomials are given corresponds to the indices of `M`. E.g., for a matrix describing a first derivative
+```
+l_i = lagrange_x
+l_j = d_lagrange_dx
+P = finite_element_matrix(l_i,l_j,0,coordinate)
+```
+the number `P[i,j]` corresponds to the matrix element
+```math
+P_{ij} = \int^{x_u}_{x_l} \Phi_i(x) \frac{d \Phi_j}{d x} d x.
+```
 For a stiffness matrix describing a second derivative, we could
 use the following.
 ```
@@ -108,13 +122,36 @@ To specify an arbitrary kernel function $`\rho(x)`$, the following function call
 ```
 coordinate = ElementCoordinates(x,scale,shift)
 M = finite_element_matrix(lagrange_x,lagrange_x,coordinate,
-            kernel_function=(v -> rho(v)),additional_quadrature_points=n)
+            kernel_function=(v -> rho(v)),additional_quadrature_points=n,
+            quadrature_increment=m,
+            max_iterations=q, atol=atol, rtol=rtol)
+```
+where we specify some function `rho(v)` to be the kernel, and we specify that we `n` additional quadrature
+points should be used in addition to the number assumed from the size of the reference grid `x`. If `q` is greater than zero, for each iteration we proceed to calculate the matrix again with a further additional `m` quadrature points, until the matrix is calculated to the specified absolute tolerance `atol` and relative tolerance `rtol`, or we call `error()` with an appropriate message for the user. Note that `rho(v)` should be specified in terms of the physical coordinate including the scale and shift factors `v = scale * x + shift`, not in terms of the reference coordinate `x`.
+
+# Matrices in two coordinates (2D)
+
+It is possible to write down weak forms where the multi-dimensional integrals do not separate into a product of one-dimensional integrals, even for a tensor-product basis. This occurs where the kernel is not a separable function, e.g., for the matrix
+```math
+K_{ikjn} = \int^{x_{2u}}_{x_{2l}}\int^{x_{1u}}_{x_{1l}} P_i(x_1) R_k(x_2) Q_j(x_1) S_n(x_2)\rho(x_1,x_2) d x_1 d x_2,
 ```
-where we specify some function `rho(v)` to be the kernel, and we specify than `n` additional quadrature
-points should be used in addition to the number assumed from the size of the reference grid `x`.
-Note that `rho(v)` should be specified in terms of the physical coordinate including the scale and shift
-factors `v = scale * x + shift`, not in terms of the reference coordinate `x`.
+where $`P_i(x_1)`$, $`R_i(x_2)`$, $`Q_i(x_1)`$, $`S_i(x_2)`$ are from the set of Lagrange polynomials basis
+functions (and their first derivatives) in $`x_1`$ and $`x_2`$, respectively.
 
+The syntax that we use to create matrices in two coordinates is as follows
+```
+l_i = lagrange_x
+l_j = lagrange_x
+l_k = d_lagrange_dx
+l_n = d_lagrange_dx
+K_2D = finite_element_matrix(l_i,l_j,coordinate_x1,
+                            l_k,l_n,coordinate_x2;
+                            kernel_function=((x1,x2)->rho(x1,x2)),
+                            max_iterations=10, quadrature_increment=10,
+                            atol=1.0e-13, rtol=1.0e-13)
+```
+where we have two coordinates, `x1` (with polynomials `l_i` and `l_j`) and `x2` (with polynomials `l_k` and `l_n`). Note that the `@enum` variables used to select the Lagrange polynomials are the same between the two coordinates.
+The matrix created above is indexed in the order `K_2D[i,k,j,n]`, indexing over the first polynomial `l_i` of `x1`, then the first polynomial `l_k` of `x2`, then the second polynomial `l_j` of `x1` with the final index representing the polynomial `l_n` of `x2`. The same adaptive quadrature keyword arguments are used as in the 1D case.
 # Examples
 
 There are several examples of taking first and second derivatives
diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl
index 39234ac..4059c02 100644
--- a/src/FiniteElementMatrices.jl
+++ b/src/FiniteElementMatrices.jl
@@ -4,7 +4,8 @@ module FiniteElementMatrices
 
 using LagrangePolynomials: lagrange_poly,
                            lagrange_poly_derivative,
-                           LagrangePolyData
+                           LagrangePolyData,
+                           LagrangePolyDataBase
 using FastGaussQuadrature: gausslegendre
 
 export lagrange_x,
@@ -49,37 +50,80 @@ struct ElementCoordinates
     end
 end
 
-function select_lagrange_function(fn_type::LagrangeFunctionType)
+function select_lagrange_function(fn_type::LagrangeFunctionType,scale::Float64)
+    # including the normalisation factors due to v = scale * x + shift
     if fn_type == lagrange_x
-        func = lagrange_poly
+        func = ((lpoly_data::LagrangePolyDataBase,z::Float64) -> lagrange_poly(lpoly_data,z))
     elseif fn_type == d_lagrange_dx
-        func = lagrange_poly_derivative
+        # 1/scale due to scale in x in d l(z(x)) /d x
+        func = ((lpoly_data::LagrangePolyDataBase,z::Float64) -> (1.0/scale)*lagrange_poly_derivative(lpoly_data,z))
     end
     return func
 end
 
-function select_lagrange_prefactor(fn_type::LagrangeFunctionType,
-                                    scale::Float64)
-    if fn_type == lagrange_x
-        prefactor = 1.0
-    elseif fn_type == d_lagrange_dx
-        prefactor = 1.0/scale
+function check_power(power::Int64, varname::String)
+    if power < 0
+        error("Require power > 0, whereas $varname=$power was passed")
     end
-    return prefactor
+    return nothing
+end
+# Method for p-adaptive quadrature in 1D and 2D.
+# Note that the the rate at which the number of quadrature
+# points is increased is the same for each coordinate.
+function finite_element_matrix(args...;
+    atol=1.0e-12::Float64, rtol=1.0e-12::Float64, max_iterations=1::Int64,
+    quadrature_increment=5::Int64, verbose=false::Bool, kwargs...)
+    matrix = _finite_element_matrix(args...;
+                adaptive_quadrature_points=0,
+                kwargs...)
+    iteration = 0
+    maxvalue = maximum(abs.(matrix))
+    maxerr = maxvalue
+    tolerance_exceeded = max_iterations > 0
+    while iteration < max_iterations && tolerance_exceeded
+        # increment the iteration
+        iteration += 1
+        matrixnew = _finite_element_matrix(args...;
+                        adaptive_quadrature_points=quadrature_increment*iteration,
+                        kwargs...)
+        maxvalue = maximum(abs.(matrixnew))
+        maxerr = maximum(abs.(matrixnew .- matrix))
+        tolerance_exceeded = maxerr > atol + rtol*maxvalue
+        # update the matrix
+        matrix .= matrixnew
+    end
+    if verbose && !(tolerance_exceeded)
+        if iteration > 0
+            println("finite_element_matrix() converged after $iteration iterations \n with maxerr = $maxerr, adaptive_quadrature_points=$(quadrature_increment*iteration) \n")
+        else
+            println("finite_element_matrix() returned a matrix without p-adaptive quadrature refinement.")
+        end
+    end
+    if tolerance_exceeded
+        error("Unable to achieve specified error bounds: \n
+              maxerr < atol + rtol*maxvalue \n
+              after max_iterations = $max_iterations with atol=$atol and rtol=$rtol we obtain \n
+              maxerr = $maxerr > $atol + $rtol*$maxvalue")
+    end
+    return matrix
 end
 
+# Method for computing polynomial integrals exactly
+# for matrices with two Lagrange polymomials in 1D
 function finite_element_matrix(
     fn1_type::LagrangeFunctionType,
     fn2_type::LagrangeFunctionType,
     power::Int64,
     coordinate::ElementCoordinates
     )
-    return finite_element_matrix(fn1_type, fn2_type, coordinate;
+    check_power(power, "power")
+    return _finite_element_matrix(fn1_type, fn2_type, coordinate;
                 kernel_function=((v -> v^power)),
                 additional_quadrature_points=power)
 end
 
-function finite_element_matrix(
+# The basic method for matrices with two Lagrange polymomials in 1D
+function _finite_element_matrix(
     fn1_type::LagrangeFunctionType,
     fn2_type::LagrangeFunctionType,
     coordinate::ElementCoordinates;
@@ -87,6 +131,8 @@ function finite_element_matrix(
     # rather than of the reference coordinate z on [-1,1]
     kernel_function=((v -> 1.0))::Function,
     additional_quadrature_points=0::Int64,
+    # argument to permit single implementation of
+    adaptive_quadrature_points=0::Int64
     )
     lpoly_data = coordinate.lpoly_data
     ngrid = length(coordinate.lpoly_data.x_nodes)
@@ -95,13 +141,11 @@ function finite_element_matrix(
     # the finite element array to be returned
     matrix = zeros(Float64,ngrid,ngrid)
     # the function objects for the required polynomials
-    lagrange1 = select_lagrange_function(fn1_type)
-    lagrange2 = select_lagrange_function(fn2_type)
-    # the normalisation factors due to v = scale * x + shift
-    prefactor1 = select_lagrange_prefactor(fn1_type,scale)
-    prefactor2 = select_lagrange_prefactor(fn2_type,scale)
-    # nquad chosen for exact results for all power, ngrid
-    nquad = ngrid + additional_quadrature_points
+    lagrange1 = select_lagrange_function(fn1_type,scale)
+    lagrange2 = select_lagrange_function(fn2_type,scale)
+    # nquad chosen for exact results for default inputs
+    # with kernel = 1.0 and zero additional quadrature points
+    nquad = ngrid + additional_quadrature_points + adaptive_quadrature_points
     zz, wz = gausslegendre(nquad)
     # compute integral
     # int P_i(z) Q_j(z) poly(z) s d z
@@ -116,14 +160,16 @@ function finite_element_matrix(
                 zzl = zz[l]
                 funcz = kernel_function(scale*zzl+shift)::Float64
                 matrix[i,j] += (scale*wz[l]*funcz*
-                           prefactor1*lagrange1(ith_lpoly_data,zzl)*
-                           prefactor2*lagrange2(jth_lpoly_data,zzl))
+                           lagrange1(ith_lpoly_data,zzl)*
+                           lagrange2(jth_lpoly_data,zzl))
             end
         end
     end
     return matrix
 end
 
+# Method for computing polynomial integrals exactly
+# for matrices with three Lagrange polymomials in 1D
 function finite_element_matrix(
     fn1_type::LagrangeFunctionType,
     fn2_type::LagrangeFunctionType,
@@ -131,11 +177,14 @@ function finite_element_matrix(
     power::Int64,
     coordinate::ElementCoordinates
     )
-    return finite_element_matrix(fn1_type, fn2_type, fn3_type,
+    check_power(power, "power")
+    return _finite_element_matrix(fn1_type, fn2_type, fn3_type,
                 coordinate; kernel_function=((v -> v^power)),
                 additional_quadrature_points=power)
 end
-function finite_element_matrix(
+
+# The basic method for matrices with three Lagrange polymomials in 2D
+function _finite_element_matrix(
     fn1_type::LagrangeFunctionType,
     fn2_type::LagrangeFunctionType,
     fn3_type::LagrangeFunctionType,
@@ -144,6 +193,7 @@ function finite_element_matrix(
     # rather than of the reference coordinate z on [-1,1]
     kernel_function=((v -> 1.0))::Function,
     additional_quadrature_points=0::Int64,
+    adaptive_quadrature_points=0::Int64
     )
     lpoly_data = coordinate.lpoly_data
     ngrid = length(coordinate.lpoly_data.x_nodes)
@@ -152,15 +202,12 @@ function finite_element_matrix(
     # the finite element array to be returned
     matrix = zeros(Float64,ngrid,ngrid,ngrid)
     # the function objects for the required polynomials
-    lagrange1 = select_lagrange_function(fn1_type)
-    lagrange2 = select_lagrange_function(fn2_type)
-    lagrange3 = select_lagrange_function(fn3_type)
-    # the normalisation factors due to v = scale * x + shift
-    prefactor1 = select_lagrange_prefactor(fn1_type,scale)
-    prefactor2 = select_lagrange_prefactor(fn2_type,scale)
-    prefactor3 = select_lagrange_prefactor(fn3_type,scale)
-    # nquad chosen for exact results for all power, ngrid
-    nquad = 2*ngrid + additional_quadrature_points
+    lagrange1 = select_lagrange_function(fn1_type,scale)
+    lagrange2 = select_lagrange_function(fn2_type,scale)
+    lagrange3 = select_lagrange_function(fn3_type,scale)
+    # nquad chosen for exact results for default inputs
+    # with kernel = 1.0 and zero additional quadrature points
+    nquad = 2*ngrid + additional_quadrature_points + adaptive_quadrature_points
     zz, wz = gausslegendre(nquad)
     # compute integral
     # int P_i(z) Q_j(z) S_k(z) poly(z) d z
@@ -175,11 +222,110 @@ function finite_element_matrix(
                 ith_lpoly_data = lpoly_data.lpoly_data[i]
                 for l in 1:nquad
                     zzl = zz[l]
-                    funcz = kernel_function(scale*zzl+shift)
+                    funcz = kernel_function(scale*zzl+shift)::Float64
                     matrix[i,j,k] += (scale*wz[l]*funcz*
-                            prefactor1*lagrange1(ith_lpoly_data,zzl)*
-                            prefactor2*lagrange2(jth_lpoly_data,zzl)*
-                            prefactor3*lagrange3(kth_lpoly_data,zzl))
+                            lagrange1(ith_lpoly_data,zzl)*
+                            lagrange2(jth_lpoly_data,zzl)*
+                            lagrange3(kth_lpoly_data,zzl))
+                end
+            end
+        end
+    end
+    return matrix
+end
+
+# Method for computing polynomial integrals exactly
+# for matrices with two 2D basis functions
+# (two Lagrange polymomials per coordinate) in 2D
+function finite_element_matrix(
+    fn1_x1_type::LagrangeFunctionType,
+    fn2_x1_type::LagrangeFunctionType,
+    power_x1::Int64,
+    coordinate_x1::ElementCoordinates,
+    fn1_x2_type::LagrangeFunctionType,
+    fn2_x2_type::LagrangeFunctionType,
+    power_x2::Int64,
+    coordinate_x2::ElementCoordinates
+    )
+    check_power(power_x1, "power_x1")
+    check_power(power_x2, "power_x2")
+    return _finite_element_matrix(fn1_x1_type, fn2_x1_type, coordinate_x1,
+                fn1_x2_type, fn2_x2_type, coordinate_x2;
+                kernel_function=((v1,v2) -> (v1^power_x1)*(v2^power_x2)),
+                additional_quadrature_points_x1=power_x1,
+                additional_quadrature_points_x2=power_x2)
+end
+
+# Basic method for computing matrices with two 2D basis functions
+# (two Lagrange polymomials per coordinate) in 2D
+function _finite_element_matrix(
+    fn1_x1_type::LagrangeFunctionType,
+    fn2_x1_type::LagrangeFunctionType,
+    coordinate_x1::ElementCoordinates,
+    fn1_x2_type::LagrangeFunctionType,
+    fn2_x2_type::LagrangeFunctionType,
+    coordinate_x2::ElementCoordinates;
+    # function of the "physical" coord v = s z + c
+    # rather than of the reference coordinate z on [-1,1]
+    kernel_function=((v1,v2) -> 1.0)::Function,
+    additional_quadrature_points_x1=0::Int64,
+    additional_quadrature_points_x2=0::Int64,
+    adaptive_quadrature_points=0::Int64
+    )
+    # coordinate x1 data
+    lpoly_data_x1 = coordinate_x1.lpoly_data
+    ngrid_x1 = length(coordinate_x1.lpoly_data.x_nodes)
+    scale_x1 = coordinate_x1.scale
+    shift_x1 = coordinate_x1.shift
+    # coordinate x2 data
+    lpoly_data_x2 = coordinate_x2.lpoly_data
+    ngrid_x2 = length(coordinate_x2.lpoly_data.x_nodes)
+    scale_x2 = coordinate_x2.scale
+    shift_x2 = coordinate_x2.shift
+    # the finite element array to be returned
+    matrix = zeros(Float64,ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
+    # the function objects for the required polynomials
+    lagrange11 = select_lagrange_function(fn1_x1_type,scale_x1)
+    lagrange12 = select_lagrange_function(fn1_x2_type,scale_x2)
+    lagrange21 = select_lagrange_function(fn2_x1_type,scale_x1)
+    lagrange22 = select_lagrange_function(fn2_x2_type,scale_x2)
+    # the normalisation factors due to v = scale * x + shift
+    # nquad chosen for exact results for default inputs
+    # with kernel = 1.0 and zero additional quadrature points
+    nquad_x1 = ngrid_x1 + additional_quadrature_points_x1 + adaptive_quadrature_points
+    zz_x1, wz_x1 = gausslegendre(nquad_x1)
+    nquad_x2 = ngrid_x2 + additional_quadrature_points_x2 + adaptive_quadrature_points
+    zz_x2, wz_x2 = gausslegendre(nquad_x2)
+    # compute integral
+    # \int \int \left(P1_i(z_1) Q1_j(z_1) P2_i(z_1) Q2_j(z_2)
+    # kernel(s_1 z_1 + c_1, s_2 z_2 + c_2) \right) s_1 s_2 d z_1 d z_2
+    # s = scale, c = shift and
+    # P_i(z), Q_i(z) in [l_i(z), (1/s) d l_i(z) / d z]
+    for jx2 in 1:ngrid_x2
+        jx2_lpoly_data = lpoly_data_x2.lpoly_data[jx2]
+        for jx1 in 1:ngrid_x1
+            jx1_lpoly_data = lpoly_data_x1.lpoly_data[jx1]
+            for ix2 in 1:ngrid_x2
+                ix2_lpoly_data = lpoly_data_x2.lpoly_data[ix2]
+                for ix1 in 1:ngrid_x1
+                    ix1_lpoly_data = lpoly_data_x1.lpoly_data[ix1]
+                    for lx2 in 1:nquad_x2
+                        zzl2 = zz_x2[lx2]
+                        wgt2 = scale_x2*wz_x2[lx2]
+                        lagrange_factor2 = (lagrange12(ix2_lpoly_data,zzl2)*
+                                               lagrange22(jx2_lpoly_data,zzl2))
+                        for lx1 in 1:nquad_x1
+                            zzl1 = zz_x1[lx1]
+                            wgt1 = scale_x1*wz_x1[lx1]
+                            lagrange_factor1 = (lagrange11(ix1_lpoly_data,zzl1)*
+                                               lagrange21(jx1_lpoly_data,zzl1))
+                            kernel_x1_x2 = kernel_function(scale_x1*zzl1+shift_x1,
+                                                        scale_x2*zzl2+shift_x2)::Float64
+                            matrix[ix1,ix2,jx1,jx2] += (wgt1*wgt2*kernel_x1_x2*
+                                                        lagrange_factor1*
+                                                        lagrange_factor2)
+                        end
+                    end
                 end
             end
         end
@@ -187,4 +333,5 @@ function finite_element_matrix(
     return matrix
 end
 
+
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 1706ed3..9299d0a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -379,6 +379,99 @@ function test_nonlinear_operator_nonpolynomial_kernel(;nodes::node_type=GLL,
     #println(maxerr)
 end
 
+function test_2D_linear_operators(;nodes::node_type=GLL,
+            ngrid_x1=5::Int64,ngrid_x2=4::Int64,
+            atol=1.0e-13::Float64)
+    # x1
+    # endpoints of grid
+    y1 = 1.2
+    y0 = -1.3
+    x1 = reference_nodes(nodes,ngrid_x1)
+    scale_x1 = 0.5*(y1-y0)
+    shift_x1 = 0.5*(y0+y1)
+    coordinate_x1 = ElementCoordinates(x1,scale_x1,shift_x1)
+    # x2
+    # endpoints of grid
+    y1 = 0.9
+    y0 = -0.78
+    x2 = reference_nodes(nodes,ngrid_x2)
+    scale_x2 = 0.5*(y1-y0)
+    shift_x2 = 0.5*(y0+y1)
+    coordinate_x2 = ElementCoordinates(x2,scale_x2,shift_x2)
+
+    M_x1 = zeros(ngrid_x1,ngrid_x1)
+    M_x2 = zeros(ngrid_x2,ngrid_x2)
+    M_2D = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
+    M_2D_err = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
+    M_2D_exact = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
+    # test that a 2D matrix with polynomial kernel can be
+    # calculated with an identical result to the
+    # product of 1D matrices with the same kernels
+    for power1 in 0:3
+        for power2 in 0:3
+            for l1A in (lagrange_x, d_lagrange_dx)
+                for l1B in (lagrange_x, d_lagrange_dx)
+                    for l2A in (lagrange_x, d_lagrange_dx)
+                        for l2B in (lagrange_x, d_lagrange_dx)
+                            #1D matrices
+                            M_x1 .= finite_element_matrix(l1A,l1B,power1,coordinate_x1)
+                            M_x2 .= finite_element_matrix(l2A,l2B,power2,coordinate_x2)
+                            #2D matrices
+                            M_2D .= finite_element_matrix(l1A,l1B,power1,coordinate_x1,
+                                                        l2A,l2B,power2,coordinate_x2)
+                            for i2 in 1:ngrid_x2
+                                for i1 in 1:ngrid_x1
+                                    for j2 in 1:ngrid_x2
+                                        for j1 in 1:ngrid_x1
+                                            M_2D_exact[j1,j2,i1,i2] = M_x1[j1,i1]*M_x2[j2,i2]
+                                        end
+                                    end
+                                end
+                            end
+                            @. M_2D_err = abs(M_2D-M_2D_exact)
+                            @test maximum(M_2D_err) < atol
+                        end
+                    end
+                end
+            end
+        end
+    end
+    # tests 2D matrices against products of 1D matrices for a simple separable Jacobian
+    for l1A in (lagrange_x, d_lagrange_dx)
+        for l1B in (lagrange_x, d_lagrange_dx)
+            for l2A in (lagrange_x, d_lagrange_dx)
+                for l2B in (lagrange_x, d_lagrange_dx)
+                    #1D matrices
+                    M_x1 .= finite_element_matrix(l1A,l1B,coordinate_x1;
+                                kernel_function=(x1 -> exp(-x1^2)), quadrature_increment=10,
+                                max_iterations=10, atol=1.0e-13, rtol=1.0e-13)
+                    M_x2 .= finite_element_matrix(l2A,l2B,coordinate_x2;
+                                kernel_function=(x2 -> exp(-x2^2)), quadrature_increment=10,
+                                max_iterations=10, atol=1.0e-13, rtol=1.0e-13)
+                    #2D matrices
+                    M_2D .= finite_element_matrix(l1A,l1B,coordinate_x1,
+                                                l2A,l2B,coordinate_x2;
+                                                kernel_function=((x1,x2)->exp(-x1^2 - x2^2)),
+                                                max_iterations=10, quadrature_increment=10,
+                                                atol=1.0e-13, rtol=1.0e-13)
+                    for i2 in 1:ngrid_x2
+                        for i1 in 1:ngrid_x1
+                            for j2 in 1:ngrid_x2
+                                for j1 in 1:ngrid_x1
+                                    M_2D_exact[j1,j2,i1,i2] = M_x1[j1,i1]*M_x2[j2,i2]
+                                end
+                            end
+                        end
+                    end
+                    @. M_2D_err = abs(M_2D-M_2D_exact)
+                    @test maximum(M_2D_err) < atol
+                end
+            end
+        end
+    end
+    return nothing
+end
+
 function runtests()
     @testset "FiniteElementMatrices" begin
         println("test FiniteElementMatrices")
@@ -414,7 +507,6 @@ function runtests()
             end
         end
         println("test second derivative:")
-        nodes_list = [GLL, GLR, GLe]
         for nodes in nodes_list
             n = 20
             println("    -- test: $(nodes) ngrid=$(n) trig")
@@ -424,7 +516,6 @@ function runtests()
         end
         
         println("test nonlinear operators:")
-        nodes_list = [GLL, GLR, GLe]
         for nodes in nodes_list
             n = 25
             println("    -- test: $(nodes) ngrid=$(n) trig")
@@ -434,13 +525,22 @@ function runtests()
         end
 
         println("test nonpolynomial kernels:")
-        nodes_list = [GLL, GLR, GLe]
         for nodes in nodes_list
             n = 20
             println("    -- test: $(nodes) ngrid=$(n) trig")
             test_nonpolynomial_kernel(; nodes=nodes, ngrid=n)
             test_nonlinear_operator_nonpolynomial_kernel(; nodes=nodes, ngrid=n)
         end
+
+        println("test 2D operators:")
+        for nodes in nodes_list
+            println("    -- test: $(nodes) ngrid_x1=5 ngrid_x2=4 poly")
+            test_2D_linear_operators(;nodes=nodes,
+                ngrid_x1=5, ngrid_x2=4, atol=1.0e-13)
+            println("    -- test: $(nodes) ngrid_x1=4 ngrid_x2=5 poly")
+            test_2D_linear_operators(;nodes=nodes,
+                ngrid_x1=4, ngrid_x2=5, atol=1.0e-13)
+        end
     end
 end