diff --git a/moment_kinetics/src/calculus.jl b/moment_kinetics/src/calculus.jl
index adff79e30..ee2910dac 100644
--- a/moment_kinetics/src/calculus.jl
+++ b/moment_kinetics/src/calculus.jl
@@ -2,9 +2,14 @@
 """
 module calculus
 
+# Import moment_kinetics so that we can refer to it in docstrings
+import moment_kinetics
+
 export derivative!, second_derivative!, laplacian_derivative!
+export elementwise_indefinite_integration!
 export reconcile_element_boundaries_MPI!
 export integral
+export indefinite_integral!
 
 using ..moment_kinetics_structs: discretization_info, null_spatial_dimension_info,
                                  null_velocity_dimension_info, weak_discretization_info
@@ -30,6 +35,55 @@ Result is stored in coord.scratch_2d.
 """
 function elementwise_derivative! end
 
+"""
+"""
+function elementwise_indefinite_integration! end
+
+"""
+    indefinite_integral!(pf, f, coord, spectral)
+
+Indefinite line integral. 
+
+This function is designed to work on local-in-memory data only,
+with distributed-memory MPI not implemented here.
+A function which integrates along a line which is distributed
+in memory exists in [`moment_kinetics.em_fields`](@ref) as 
+`calculate_phi_from_Epar!()`. The distributed-memory functionality could
+be ported to a generic function, similiar to how the derivative!
+functions are generalised in [`moment_kinetics.derivatives`](@ref).
+"""
+function indefinite_integral!(pf, f, coord, spectral::discretization_info)
+    # get the indefinite integral at each grid point within each element and store in
+    # coord.scratch_2d
+    elementwise_indefinite_integration!(coord, f, spectral)
+    # map the integral from the elemental grid to the full grid;
+    # taking care to match integral constants at the boundaries
+    # we assume that the lower limit of the indefinite integral
+    # is the lowest value in coord.grid
+    indefinite_integral_elements_to_full_grid!(pf, coord)
+end
+
+"""
+"""
+function indefinite_integral_elements_to_full_grid!(pf, coord)
+    pf2d = coord.scratch_2d
+    nelement_local = coord.nelement_local
+    ngrid = coord.ngrid
+    igrid_full = coord.igrid_full   
+    j = 1 # the first element
+    for i in 1:ngrid
+        pf[i] = pf2d[i,j]
+    end
+    for j in 2:nelement_local
+        ilast = igrid_full[ngrid,j-1] # the grid index of the last point in the previous element
+        for k in 2:coord.ngrid
+            i = coord.igrid_full[k,j] # the index in the full grid
+            pf[i] = pf2d[k,j] + pf[ilast]
+        end
+    end
+    return nothing
+end
+
 """
     derivative!(df, f, coord, adv_fac, spectral)
 
diff --git a/moment_kinetics/src/chebyshev.jl b/moment_kinetics/src/chebyshev.jl
index 716ca6599..479cddc98 100644
--- a/moment_kinetics/src/chebyshev.jl
+++ b/moment_kinetics/src/chebyshev.jl
@@ -17,10 +17,11 @@ using ..type_definitions: mk_float, mk_int
 using ..array_allocation: allocate_float, allocate_complex
 using ..clenshaw_curtis: clenshawcurtisweights
 import ..calculus: elementwise_derivative!
+import ..calculus: elementwise_indefinite_integration!
 using ..communication
 import ..interpolation: single_element_interpolate!
 using ..moment_kinetics_structs: discretization_info
-
+using ..gauss_legendre: integration_matrix!
 """
 Chebyshev pseudospectral discretization
 """
@@ -43,6 +44,8 @@ struct chebyshev_base_info{TForward <: FFTW.cFFTWPlan, TBackward <: AbstractFFTs
     Dmat::Array{mk_float,2}
     # elementwise differentiation vector (ngrid) for the point x = -1
     D0::Array{mk_float,1}
+    # elementwise antidifferentiation matrix
+    Amat::Array{mk_float,2}
 end
 
 struct chebyshev_info{TForward <: FFTW.cFFTWPlan, TBackward <: AbstractFFTs.ScaledPlan} <: discretization_info
@@ -141,9 +144,12 @@ function setup_chebyshev_pseudospectral_lobatto(coord, fftw_flags)
     cheb_derivative_matrix_elementwise!(Dmat,coord.ngrid)
     D0 = allocate_float(coord.ngrid)
     D0 .= Dmat[1,:]
+    Amat = allocate_float(coord.ngrid, coord.ngrid)
+    x = chebyshevpoints(coord.ngrid)
+    integration_matrix!(Amat,x,coord.ngrid)
     # return a structure containing the information needed to carry out
     # a 1D Chebyshev transform
-    return chebyshev_base_info(fext, fcheby, dcheby, forward_transform, backward_transform, Dmat, D0)
+    return chebyshev_base_info(fext, fcheby, dcheby, forward_transform, backward_transform, Dmat, D0, Amat)
 end
 
 function setup_chebyshev_pseudospectral_radau(coord, fftw_flags)
@@ -164,9 +170,12 @@ function setup_chebyshev_pseudospectral_radau(coord, fftw_flags)
         cheb_derivative_matrix_elementwise_radau_by_FFT!(Dmat, coord, fcheby, dcheby, fext, forward_transform)
         D0 = allocate_float(coord.ngrid)
         cheb_lower_endpoint_derivative_vector_elementwise_radau_by_FFT!(D0, coord, fcheby, dcheby, fext, forward_transform)
+        Amat = allocate_float(coord.ngrid, coord.ngrid)
+        x = chebyshev_radau_points(coord.ngrid)
+        integration_matrix!(Amat,x,coord.ngrid)
         # return a structure containing the information needed to carry out
         # a 1D Chebyshev transform
-        return chebyshev_base_info(fext, fcheby, dcheby, forward_transform, backward_transform, Dmat, D0)
+        return chebyshev_base_info(fext, fcheby, dcheby, forward_transform, backward_transform, Dmat, D0, Amat)
 end
 
 """
@@ -190,11 +199,10 @@ imax -- the array of maximum indices of each element on the extended grid.
 """
 function scaled_chebyshev_grid(ngrid, nelement_local, n,
 			element_scale, element_shift, imin, imax)
-    # initialize chebyshev grid defined on [1,-1]
+    # initialize chebyshev grid defined on [-1,1]
     # with n grid points chosen to facilitate
     # the fast Chebyshev transform (aka the discrete cosine transform)
     # needed to obtain Chebyshev spectral coefficients
-    # this grid goes from +1 to -1
     chebyshev_grid = chebyshevpoints(ngrid)
     # create array for the full grid
     grid = allocate_float(n)
@@ -206,9 +214,9 @@ function scaled_chebyshev_grid(ngrid, nelement_local, n,
     @inbounds for j ∈ 1:nelement_local
         scale_factor = element_scale[j]
         shift = element_shift[j]
-        # reverse the order of the original chebyshev_grid (ran from [1,-1])
+        # take the original chebyshev_grid (ran from [-1,1])
         # and apply the scale factor and shift
-        grid[imin[j]:imax[j]] .= (reverse(chebyshev_grid)[k:ngrid] * scale_factor) .+ shift
+        grid[imin[j]:imax[j]] .= (chebyshev_grid[k:ngrid] * scale_factor) .+ shift
         # after first element, increase minimum index for chebyshev_grid to 2
         # to avoid double-counting boundary element
         k = 2
@@ -219,11 +227,10 @@ end
 
 function scaled_chebyshev_radau_grid(ngrid, nelement_local, n,
 			element_scale, element_shift, imin, imax, irank)
-    # initialize chebyshev grid defined on [1,-1]
+    # initialize chebyshev grid defined on [-1,1]
     # with n grid points chosen to facilitate
     # the fast Chebyshev transform (aka the discrete cosine transform)
     # needed to obtain Chebyshev spectral coefficients
-    # this grid goes from +1 to -1
     chebyshev_grid = chebyshevpoints(ngrid)
     chebyshev_radau_grid = chebyshev_radau_points(ngrid)
     # create array for the full grid
@@ -242,9 +249,9 @@ function scaled_chebyshev_radau_grid(ngrid, nelement_local, n,
         @inbounds for j ∈ 2:nelement_local
             scale_factor = element_scale[j]
             shift = element_shift[j]
-            # reverse the order of the original chebyshev_grid (ran from [1,-1])
+            # take the original chebyshev_grid (ran from [-1,1])
             # and apply the scale factor and shift
-            grid[imin[j]:imax[j]] .= (reverse(chebyshev_grid)[k:ngrid] * scale_factor) .+ shift
+            grid[imin[j]:imax[j]] .= (chebyshev_grid[k:ngrid] * scale_factor) .+ shift
         end
         wgts = clenshaw_curtis_radau_weights(ngrid, nelement_local, n, imin, imax, element_scale)
     else
@@ -255,9 +262,9 @@ function scaled_chebyshev_radau_grid(ngrid, nelement_local, n,
         @inbounds for j ∈ 1:nelement_local
             scale_factor = element_scale[j]
             shift = element_shift[j]
-            # reverse the order of the original chebyshev_grid (ran from [1,-1])
+            # take the original chebyshev_grid (ran from [-1,1])
             # and apply the scale factor and shift
-            grid[imin[j]:imax[j]] .= (reverse(chebyshev_grid)[k:ngrid] * scale_factor) .+ shift
+            grid[imin[j]:imax[j]] .= (chebyshev_grid[k:ngrid] * scale_factor) .+ shift
             # after first element, increase minimum index for chebyshev_grid to 2
             # to avoid double-counting boundary element
             k = 2
@@ -621,6 +628,7 @@ end
 
 """
 returns the Chebyshev-Gauss-Lobatto grid points on an n point grid
+in the range [-1,1]
 """
 function chebyshevpoints(n)
     grid = allocate_float(n)
@@ -631,9 +639,14 @@ function chebyshevpoints(n)
             grid[j] = cospi((j-1)*nfac)
         end
     end
-    return grid
+    # return grid on z in [-1,1]
+    return reverse(grid)
 end
 
+"""
+returns the Chebyshev-Gauss-Radau grid points on an n point grid
+in the range (-1,1]
+"""
 function chebyshev_radau_points(n)
     grid = allocate_float(n)
     nfac = 1.0/(n-0.5)
@@ -979,4 +992,46 @@ https://people.maths.ox.ac.uk/trefethen/pdetext.html
         D[1] = -sum(D[:])
     end
 
+"""
+A function that takes the indefinite integral in each element of `coord.grid`,
+leaving the result (element-wise) in `coord.scratch_2d`.
+"""
+function elementwise_indefinite_integration!(coord, ff, chebyshev::chebyshev_info)
+    # the primitive of f
+    pf = coord.scratch_2d
+    # define local variable nelement for convenience
+    nelement = coord.nelement_local
+    # check array bounds
+    @boundscheck nelement == size(pf,2) && coord.ngrid == size(pf,1) || throw(BoundsError(pf))
+    
+    # variable k will be used to avoid double counting of overlapping point
+    k = 0
+    j = 1 # the first element
+    imin = coord.imin[j]-k
+    # imax is the maximum index on the full grid for this (jth) element
+    imax = coord.imax[j]        
+    if coord.radau_first_element && coord.irank == 0 # differentiate this element with the Radau scheme
+        @views mul!(pf[:,j],chebyshev.radau.Amat[:,:],ff[imin:imax])
+    else #differentiate using the Lobatto scheme
+        @views mul!(pf[:,j],chebyshev.lobatto.Amat[:,:],ff[imin:imax])
+    end
+    # transform back to the physical coordinate scale
+    for i in 1:coord.ngrid
+        pf[i,j] *= coord.element_scale[j]
+    end
+    # calculate the derivative on each element
+    @inbounds for j ∈ 2:nelement
+        k = 1 
+        imin = coord.imin[j]-k
+        # imax is the maximum index on the full grid for this (jth) element
+        imax = coord.imax[j]
+        @views mul!(pf[:,j],chebyshev.lobatto.Amat[:,:],ff[imin:imax])        
+        # transform back to the physical coordinate scale
+        for i in 1:coord.ngrid
+            pf[i,j] *= coord.element_scale[j]
+        end
+    end
+    return nothing
+end
+
 end
diff --git a/moment_kinetics/src/em_fields.jl b/moment_kinetics/src/em_fields.jl
index 462d23959..3a7643912 100644
--- a/moment_kinetics/src/em_fields.jl
+++ b/moment_kinetics/src/em_fields.jl
@@ -4,6 +4,8 @@ module em_fields
 
 export setup_em_fields
 export update_phi!
+# for testing
+export calculate_phi_from_Epar!
 
 using ..type_definitions: mk_float
 using ..array_allocation: allocate_shared_float
@@ -16,6 +18,7 @@ using ..velocity_moments: update_density!
 using ..derivatives: derivative_r!, derivative_z!
 using ..electron_fluid_equations: calculate_Epar_from_electron_force_balance!
 using ..gyroaverages: gyro_operators, gyroaverage_field!
+using ..calculus: indefinite_integral!
 
 using MPI
 
@@ -121,7 +124,7 @@ function update_phi!(fields, fvec, vperp, z, r, composition, collisions, moments
             composition.n_neutral_species, collisions.reactions.electron_charge_exchange_frequency,
             composition.me_over_mi, fvec.density_neutral, fvec.uz_neutral,
             fvec.electron_upar)
-        calculate_phi_from_Epar!(fields.phi, fields.Ez, r, z)
+        calculate_phi_from_Epar!(fields.phi, fields.Ez, r, z, z_spectral)
     end
     ## can calculate phi at z = L and hence phi_wall(z=L) using jpar_i at z =L if needed
     _block_synchronize()
@@ -174,34 +177,49 @@ function update_phi!(fields, fvec, vperp, z, r, composition, collisions, moments
     return nothing
 end
 
-function calculate_phi_from_Epar!(phi, Epar, r, z)
+function calculate_phi_from_Epar!(phi, Epar, r, z, z_spectral)
     # simple calculation of phi from Epar for now. Assume phi is already set at the
     # lower-ze boundary, e.g. by the kinetic electron boundary condition.
     begin_serial_region()
 
-    dz = z.cell_width
     @serial_region begin
-        # Need to broadcast the lower-z boundary value, because we only communicate
-        # delta_phi below, rather than passing the boundary values directly from block to
-        # block.
-        MPI.Bcast!(@view(phi[1,:]), z.comm; root=0)
-
+        if z.irank == 0
+            # Don't want to change the lower-z boundary value, so save it here so we can
+            # restore it at the end
+            @views @. r.scratch3 = phi[1,:]
+        end
+        # Broadcast the lower-z boundary value, so that we can communicate
+        # delta_phi = phi[end,:] - phi[1,:] below, rather than passing the boundary values directly
+        # from block to block.
+        MPI.Bcast!(r.scratch3, z.comm; root=0)
         if z.irank == z.nrank - 1
             # Don't want to change the upper-z boundary value, so save it here so we can
             # restore it at the end
             @views @. r.scratch = phi[end,:]
         end
-
+        
+        # calculate phi on each local rank, up to a constant
+        # phi[1,:] = 0.0 by convention here
         @loop_r ir begin
-            for iz ∈ 2:z.n
-                phi[iz,ir] = phi[iz-1,ir] - dz[iz-1]*Epar[iz,ir]
-            end
+            @views @. z.scratch = -Epar[:,ir]
+            @views indefinite_integral!(phi[:,ir], z.scratch, z, z_spectral)
+        end
+        
+        # Restore the constant offset from the lower boundary
+        # to all ranks so that we can use this_delta_phi below
+        @loop_r_z ir iz begin
+            @views phi[iz,ir] += r.scratch3[ir]
         end
-
         # Add contributions to integral along z from processes at smaller z-values than
         # this one.
-        this_delta_phi = r.scratch2 .= phi[end,:] .- phi[1,:]
+        this_delta_phi = r.scratch2
         for irank ∈ 0:z.nrank-2
+            # update this_delta_phi before broadcasting
+            @loop_r ir begin
+                this_delta_phi[ir] = phi[end,ir] - phi[1,ir]
+            end
+            # broadcast this_delta_phi from z.irank = irank to all processes
+            # updating this_delta_phi everywhere
             MPI.Bcast!(this_delta_phi, z.comm; root=irank)
             if z.irank > irank
                 @loop_r_z ir iz begin
@@ -212,7 +230,9 @@ function calculate_phi_from_Epar!(phi, Epar, r, z)
 
         if z.irank == z.nrank - 1
             # Restore the upper-z boundary value
-            @views @. phi[end,:] = r.scratch
+            @loop_r ir begin 
+                @views phi[end,ir] = r.scratch[ir]
+            end
         end
     end
 
diff --git a/moment_kinetics/src/finite_differences.jl b/moment_kinetics/src/finite_differences.jl
index cc2d9583b..845298f0e 100644
--- a/moment_kinetics/src/finite_differences.jl
+++ b/moment_kinetics/src/finite_differences.jl
@@ -5,6 +5,7 @@ module finite_differences
 using ..type_definitions: mk_float
 import ..calculus: elementwise_derivative!, second_derivative!,
                    derivative_elements_to_full_grid!
+import ..calculus: elementwise_indefinite_integration!
 using ..moment_kinetics_structs: discretization_info
 
 """
@@ -540,4 +541,51 @@ function second_derivative_finite_difference!(df::Array{mk_float,2}, f, del, bc,
     end
 end
 
+"""
+    elementwise_indefinite_integration!(coord, f, not_spectral::finite_difference_info)
+
+Calculate the primitive of f using second-order accurate trapezium rule; result stored
+in coord.scratch_2d.
+"""
+function elementwise_indefinite_integration!(coord, f, not_spectral::finite_difference_info)
+    return primitive_finite_difference!(coord, f, coord.finite_difference_option)
+end
+
+"""
+"""
+function primitive_finite_difference!(coord, f, finite_difference_option)
+	# space here to add different integration methods, if required
+	primitive_second_order!(coord, f)
+	return nothing
+end
+
+"""
+Integrate the input function f and return as pf
+using second-order trapezium rule.
+Do the integral on each element separately.
+"""
+function primitive_second_order!(coord, f)
+	n = length(f)
+	ngrid = coord.ngrid
+	nelement = coord.nelement_local
+	# the primitive
+	pf = coord.scratch_2d
+	# the grid
+	z = coord.grid
+	igrid_full = coord.igrid_full
+	for j in 1:nelement
+		# lower value of primitive is zero
+		# indefinite_integral_elements_to_full_grid!() in calculus.jl
+		# will calculate the correct integration constants
+		pf[1,j] = 0.0
+		# do the integral on the element
+		for i in 2:ngrid
+			# the index on the full grid
+			k = igrid_full[i,j]
+			# the summation
+			pf[i,j] = pf[i-1,j] + 0.5*(z[k] - z[k-1])*(f[k]+f[k-1])
+		end
+    end
+end
+
 end
diff --git a/moment_kinetics/src/gauss_legendre.jl b/moment_kinetics/src/gauss_legendre.jl
index bb5135ba9..8dcbfb176 100644
--- a/moment_kinetics/src/gauss_legendre.jl
+++ b/moment_kinetics/src/gauss_legendre.jl
@@ -20,6 +20,7 @@ export setup_gausslegendre_pseudospectral
 export GaussLegendre_weak_product_matrix!
 export ielement_global_func
 export get_QQ_local!
+export integration_matrix!
 
 using FastGaussQuadrature
 using LegendrePolynomials: Pl, dnPl
@@ -28,10 +29,11 @@ using SparseArrays: sparse, AbstractSparseArray
 using SparseMatricesCSR
 using ..type_definitions: mk_float, mk_int
 using ..array_allocation: allocate_float
-import ..calculus: elementwise_derivative!, mass_matrix_solve!
+import ..calculus: elementwise_derivative!, mass_matrix_solve!,
+                   elementwise_indefinite_integration!
 import ..interpolation: single_element_interpolate!,
                         fill_single_element_interpolation_matrix!
-using ..lagrange_polynomials: lagrange_poly_optimised, lagrange_poly_derivative_optimised
+using ..lagrange_polynomials: lagrange_poly_optimised, lagrange_poly_derivative_optimised, lagrange_poly
 using ..moment_kinetics_structs: weak_discretization_info
 
 
@@ -44,6 +46,8 @@ on Gauss-Legendre points in 1D
 struct gausslegendre_base_info
     # elementwise differentiation matrix (ngrid*ngrid)
     Dmat::Array{mk_float,2}
+    # elementwise integration matrix (ngrid*ngrid)
+    Amat::Array{mk_float,2}
     # local mass matrix type 0
     M0::Array{mk_float,2}
     # local mass matrix type 1
@@ -186,6 +190,8 @@ function setup_gausslegendre_pseudospectral_lobatto(coord; collision_operator_di
     w = mk_float.(w)
     Dmat = allocate_float(coord.ngrid, coord.ngrid)
     gausslobattolegendre_differentiation_matrix!(Dmat,x,coord.ngrid)
+    Amat = allocate_float(coord.ngrid, coord.ngrid)
+    integration_matrix!(Amat,x,coord.ngrid)
     
     M0 = allocate_float(coord.ngrid, coord.ngrid)
     GaussLegendre_weak_product_matrix!(M0,coord.ngrid,x,w,"M0")
@@ -248,7 +254,7 @@ function setup_gausslegendre_pseudospectral_lobatto(coord; collision_operator_di
         Y30 = allocate_float(0, 0, 0)
         Y31 = allocate_float(0, 0, 0)
     end
-    return gausslegendre_base_info(Dmat,M0,M1,M2,S0,S1,
+    return gausslegendre_base_info(Dmat,Amat,M0,M1,M2,S0,S1,
             K0,K1,K2,P0,P1,P2,D0,Y00,Y01,Y10,Y11,Y20,Y21,Y30,Y31)
 end
 
@@ -267,7 +273,9 @@ function setup_gausslegendre_pseudospectral_radau(coord; collision_operator_dim=
     # elemental differentiation matrix
     Dmat = allocate_float(coord.ngrid, coord.ngrid)
     gaussradaulegendre_differentiation_matrix!(Dmat,x,coord.ngrid)
-    
+    Amat = allocate_float(coord.ngrid, coord.ngrid)
+    integration_matrix!(Amat,xreverse,coord.ngrid)
+
     M0 = allocate_float(coord.ngrid, coord.ngrid)
     GaussLegendre_weak_product_matrix!(M0,coord.ngrid,xreverse,wreverse,"M0",radau=true)
     M1 = allocate_float(coord.ngrid, coord.ngrid)
@@ -327,9 +335,50 @@ function setup_gausslegendre_pseudospectral_radau(coord; collision_operator_dim=
         Y30 = allocate_float(0, 0, 0)
         Y31 = allocate_float(0, 0, 0)
     end
-    return gausslegendre_base_info(Dmat,M0,M1,M2,S0,S1,
+    return gausslegendre_base_info(Dmat,Amat,M0,M1,M2,S0,S1,
             K0,K1,K2,P0,P1,P2,D0,Y00,Y01,Y10,Y11,Y20,Y21,Y30,Y31)
 end 
+"""
+A function that takes the indefinite integral in each element of `coord.grid`,
+leaving the result (element-wise) in `coord.scratch_2d`.
+"""
+function elementwise_indefinite_integration!(coord, ff, gausslegendre::gausslegendre_info)
+    # the primitive of f
+    pf = coord.scratch_2d
+    # define local variable nelement for convenience
+    nelement = coord.nelement_local
+    # check array bounds
+    @boundscheck nelement == size(pf,2) && coord.ngrid == size(pf,1) || throw(BoundsError(pf))
+    
+    # variable k will be used to avoid double counting of overlapping point
+    k = 0
+    j = 1 # the first element
+    imin = coord.imin[j]-k
+    # imax is the maximum index on the full grid for this (jth) element
+    imax = coord.imax[j]        
+    if coord.radau_first_element && coord.irank == 0 # differentiate this element with the Radau scheme
+        @views mul!(pf[:,j],gausslegendre.radau.Amat[:,:],ff[imin:imax])
+    else #differentiate using the Lobatto scheme
+        @views mul!(pf[:,j],gausslegendre.lobatto.Amat[:,:],ff[imin:imax])
+    end
+    # transform back to the physical coordinate scale
+    for i in 1:coord.ngrid
+        pf[i,j] *= coord.element_scale[j]
+    end
+    # calculate the derivative on each element
+    @inbounds for j ∈ 2:nelement
+        k = 1 
+        imin = coord.imin[j]-k
+        # imax is the maximum index on the full grid for this (jth) element
+        imax = coord.imax[j]
+        @views mul!(pf[:,j],gausslegendre.lobatto.Amat[:,:],ff[imin:imax])        
+        # transform back to the physical coordinate scale
+        for i in 1:coord.ngrid
+            pf[i,j] *= coord.element_scale[j]
+        end
+    end
+    return nothing
+end
 
 """
 A function that takes the first derivative in each element of `coord.grid`,
@@ -463,6 +512,65 @@ function mass_matrix_solve!(f, b, spectral::gausslegendre_info)
     return nothing
 end
 
+"""
+Function to calculate the elemental anti-differentiation (integration) matrix.
+This function forms the primitive
+```math
+F(x) = \\int^x_{x_{\\rm min}} f(x^\\prime) d x^\\prime
+```
+of the function 
+```math
+f(x) = \\sum_j f_j l_j(x),
+```
+where \$l_j(x)\$ is the \$j^{\\rm th}\$ Lagrange polynomial on the element and \$f_j = f(x_j)\$,
+with \$x_j\$ \$j^{\\rm th}\$ collocation point on the element. We find \$F(x)\$ at the collocation
+points on the element, giving a series of integrals to evaluate:
+```math
+F(x_i) = \\int^{x_i}_{-1} f(x^\\prime) d x^\\prime = \\sum_j f_j \\int^{x_i}_{-1} l_j(x^\\prime) d x^\\prime,
+```
+where we have used that \$x_{\\rm min} = -1\$ on the elemental grid.
+Changing to a normalised coordinate \$y\$ suitable for Gaussian quadrature 
+```math
+x^\\prime = \\frac{x_i + 1}{2} y + \\frac{x_i - 1}{2}
+```
+we can write the operation in matrix form:
+```math
+F(x_i) = \\sum_{j}A_{ij}f_j,
+```
+with the matrix \$A_{ij}\$ defined by
+```math
+A_{ij} = \\left(\\frac{x_i + 1}{2}\\right) \\int^1_{-1} l_j \\left( \\frac{(x_i + 1)y + x_i - 1}{2} \\right) dy,
+```
+or in discretised form
+```math
+A_{ij} = \\left(\\frac{x_i + 1}{2}\\right) \\sum_k l_j \\left( \\frac{(x_i + 1)y_k + x_i - 1}{2} \\right) w_k,
+```
+with \$y_k\$ and \$w_k\$ Gauss-quadrature points and weights, respectively.
+"""
+function integration_matrix!(A::Array{Float64,2},x::Array{Float64,1},ngrid::Int64)
+    nquad = 2*ngrid
+    zz, wz = gausslegendre(nquad)
+    # set lower limit
+    for j in 1:ngrid
+        A[1,j] = 0.0
+    end
+    # set the remaining points
+    for i in 1:ngrid
+        for j in 1:ngrid
+            xi = x[i] # limit of ith integral
+            scale = 0.5*(xi + 1)
+            shift = 0.5*(xi - 1)
+            # calculate the sum with Gaussian quadrature
+            A[i,j] = 0.0
+            for k in 1:nquad
+                xval = scale*zz[k] + shift
+                A[i,j] += scale*lagrange_poly(j,x,xval)*wz[k]
+            end
+        end
+    end
+    return nothing
+end
+
 """
 Formula for Gauss-Legendre-Lobatto differentiation matrix taken from p196 of Chpt `The Spectral Elemtent Method' of 
 `Computational Seismology'. Heiner Igel First Edition. Published in 2017 by Oxford University Press.
diff --git a/moment_kinetics/test/calculus_tests.jl b/moment_kinetics/test/calculus_tests.jl
index 23a77e71a..4e55a7bbe 100644
--- a/moment_kinetics/test/calculus_tests.jl
+++ b/moment_kinetics/test/calculus_tests.jl
@@ -4,11 +4,18 @@ include("setup.jl")
 
 using moment_kinetics.coordinates: define_test_coordinate
 using moment_kinetics.calculus: derivative!, second_derivative!, integral
-using moment_kinetics.calculus: laplacian_derivative!
+using moment_kinetics.calculus: laplacian_derivative!, indefinite_integral!
 
 using MPI
 using Random
 using LinearAlgebra: mul!, ldiv!
+# To update one of the long list of tolerances, uncomment the line below
+# #using Printf
+# and insert a print statement, e.g.,
+# tolerance = 1.5*maximum(abs.(f .- expected_f))
+# tolerance_string = @sprintf("%.2g",tolerance)
+# println("($nelement, $ngrid, $tolerance_string),")
+# and comment out the failing test that needs updating or replacing
 
 function runtests()
     @testset "calculus" verbose=use_verbose begin
@@ -1863,6 +1870,194 @@ function runtests()
                                norm=maxabs_norm)
             end
         end
+
+        @testset "Indefinite line integration: spectral" verbose=false begin
+            @testset "$nelement $ngrid" for discretization in ("gausslegendre_pseudospectral", "chebyshev_pseudospectral"), name in ("z","vperp"), (nelement, ngrid, rtol) ∈
+                    (
+                     (1, 8, 1.e-3),
+                     (1, 9, 7.e-5),
+                     (1, 10, 5.e-5),
+                     (1, 11, 8.e-7),
+                     (1, 12, 5.e-7),
+                     (1, 13, 8.e-9),
+                     (1, 14, 5.e-9),
+                     (1, 15, 1.e-10),
+                     (1, 16, 5.e-10),
+                     (1, 17, 5.e-10),
+                     
+                     (2, 6, 2.e-4),
+                     (2, 7, 5.e-5),
+                     (2, 8, 1.e-6),
+                     (2, 9, 5.e-7),
+                     (2, 10, 5.e-9),
+                     (2, 11, 1.e-9),
+                     (2, 12, 5.e-10),
+                     (2, 13, 5.e-10),
+                     (2, 14, 5.e-10),
+                     (2, 15, 5.e-10),
+                     (2, 16, 5.e-10),
+                     (2, 17, 5.e-10),
+                     
+                     (3, 6, 5.e-5),
+                     (3, 7, 1.e-6),
+                     (3, 8, 8.e-8),
+                     (3, 9, 5.e-9),
+                     (3, 10, 5.e-10),
+                     (3, 11, 5.e-10),
+                     (3, 12, 5.e-10),
+                     (3, 13, 5.e-10),
+                     (3, 14, 5.e-10),
+                     (3, 15, 5.e-10),
+                     (3, 16, 5.e-10),
+                     (3, 17, 5.e-8),
+                     
+                     (4, 5, 6.e-5),
+                     (4, 6, 5.e-6),
+                     (4, 7, 1.e-7),
+                     (4, 8, 5.e-9),
+                     (4, 9, 1.e-9),
+                     (4, 10, 1.e-10),
+                     (4, 11, 8.e-10),
+                     (4, 12, 8.e-10),
+                     (4, 13, 8.e-10),
+                     (4, 14, 8.e-10),
+                     (4, 15, 8.e-10),
+                     (4, 16, 8.e-8),
+                     (4, 17, 8.e-8),
+                     
+                     (5, 5, 4.e-5),
+                     (5, 6, 8.e-7),
+                     (5, 7, 5.e-8),
+                     (5, 8, 1.e-9),
+                     (5, 9, 8.e-10),
+                     (5, 10, 5.e-10),
+                     (5, 11, 8.e-10),
+                     (5, 12, 4.e-10),
+                     (5, 13, 2.e-10),
+                     (5, 14, 2.e-10),
+                     (5, 15, 8.e-10),
+                     (5, 16, 8.e-10),
+                     (5, 17, 8.e-10),
+                     )
+
+                # define inputs needed for the test
+                L = 6.0
+                bc = "zero"
+                element_spacing_option = "uniform"
+                # create the coordinate struct 'x'
+                # This test runs effectively in serial, so implicitly uses
+                # `ignore_MPI=true` to avoid errors due to communicators not being fully
+                # set up.
+                x, spectral = define_test_coordinate(name; ngrid=ngrid,
+                                                     nelement=nelement, L=L,
+                                                     discretization=discretization,
+                                                     bc=bc,
+                                                     element_spacing_option=element_spacing_option,
+                                                     collision_operator_dim=false)
+                xllim = x.element_boundaries[1] # lower endpoint
+                f = @. cos(x.grid - xllim)
+                expected_pf = @. sin(x.grid - xllim)
+                # create array for the indefinite integral pf
+                pf = similar(f)
+
+                # differentiate f
+                indefinite_integral!(pf, f, x, spectral)
+                @test isapprox(pf, expected_pf, rtol=rtol, atol=1.e-10,
+                               norm=maxabs_norm)
+            end
+        end
+
+        @testset "Indefinite line integration: finite differences" verbose=false begin
+            @testset "$nelement $ngrid" for name in ("z","vperp"), (nelement, ngrid, rtol) ∈
+                    (
+                        (1, 8, 0.092),
+                        (1, 9, 0.071),
+                        (1, 10, 0.056),
+                        (1, 11, 0.045),
+                        (1, 12, 0.037),
+                        (1, 13, 0.031),
+                        (1, 14, 0.027),
+                        (1, 15, 0.023),
+                        (1, 16, 0.02),
+                        (1, 17, 0.018),
+                        (2, 6, 0.045),
+                        (2, 7, 0.031),
+                        (2, 8, 0.023),
+                        (2, 9, 0.018),
+                        (2, 10, 0.014),
+                        (2, 11, 0.011),
+                        (2, 12, 0.0093),
+                        (2, 13, 0.0078),
+                        (2, 14, 0.0067),
+                        (2, 15, 0.0057),
+                        (2, 16, 0.005),
+                        (2, 17, 0.0044),
+                        (3, 6, 0.02),
+                        (3, 7, 0.014),
+                        (3, 8, 0.01),
+                        (3, 9, 0.0078),
+                        (3, 10, 0.0062),
+                        (3, 11, 0.005),
+                        (3, 12, 0.0041),
+                        (3, 13, 0.0035),
+                        (3, 14, 0.003),
+                        (3, 15, 0.0026),
+                        (3, 16, 0.0022),
+                        (3, 17, 0.002),
+                        (4, 5, 0.018),
+                        (4, 6, 0.011),
+                        (4, 7, 0.0078),
+                        (4, 8, 0.0057),
+                        (4, 9, 0.0044),
+                        (4, 10, 0.0035),
+                        (4, 11, 0.0028),
+                        (4, 12, 0.0023),
+                        (4, 13, 0.002),
+                        (4, 14, 0.0017),
+                        (4, 15, 0.0014),
+                        (4, 16, 0.0013),
+                        (4, 17, 0.0011),
+                        (5, 5, 0.011),
+                        (5, 6, 0.0072),
+                        (5, 7, 0.005),
+                        (5, 8, 0.0037),
+                        (5, 9, 0.0028),
+                        (5, 10, 0.0022),
+                        (5, 11, 0.0018),
+                        (5, 12, 0.0015),
+                        (5, 13, 0.0013),
+                        (5, 14, 0.0011),
+                        (5, 15, 0.00092),
+                        (5, 16, 0.0008),
+                        (5, 17, 0.0007),
+                     )
+
+                # define inputs needed for the test
+                L = 6.0
+                bc = "zero"
+                element_spacing_option = "uniform"
+                # create the coordinate struct 'x'
+                # This test runs effectively in serial, so implicitly uses
+                # `ignore_MPI=true` to avoid errors due to communicators not being fully
+                # set up.
+                x, spectral = define_test_coordinate(name; ngrid=ngrid,
+                                                     nelement=nelement, L=L,
+                                                     discretization="finite_difference",
+                                                     bc=bc,
+                                                     element_spacing_option=element_spacing_option,
+                                                     collision_operator_dim=false)
+                xllim = x.element_boundaries[1] # lower endpoint
+                f = @. cos(x.grid - xllim)
+                expected_pf = @. sin(x.grid - xllim)
+                # create array for the indefinite integral pf
+                pf = similar(f)
+
+                # differentiate f
+                indefinite_integral!(pf, f, x, spectral)
+                @test isapprox(pf, expected_pf, rtol=rtol, atol=1.e-10,
+                               norm=maxabs_norm)
+            end
+        end
     end
 end
 
diff --git a/test_scripts/GaussLobattoLegendre_test.jl b/test_scripts/GaussLobattoLegendre_test.jl
index 74292aa1e..27c232726 100644
--- a/test_scripts/GaussLobattoLegendre_test.jl
+++ b/test_scripts/GaussLobattoLegendre_test.jl
@@ -11,7 +11,7 @@ import moment_kinetics
 using moment_kinetics.gauss_legendre
 using moment_kinetics.coordinates: define_coordinate
 using moment_kinetics.calculus: derivative!, second_derivative!, laplacian_derivative!
-using moment_kinetics.calculus: mass_matrix_solve!
+using moment_kinetics.calculus: mass_matrix_solve!, indefinite_integral!
 using moment_kinetics.type_definitions: OptionsDict
 
 
@@ -71,7 +71,34 @@ using moment_kinetics.type_definitions: OptionsDict
             # This test runs effectively in serial, so use `ignore_MPI=true` to avoid
             # errors due to communicators not being fully set up.
             y, y_spectral = define_coordinate(input, y_name; collision_operator_dim=true, ignore_MPI=true)
-            #print_matrix(Mmat,"Mmat",y.n,y.n)
+            print_matrix(y_spectral.lobatto.Amat,"Amat",y.ngrid,y.ngrid)
+            print_matrix(y_spectral.radau.Amat,"Amat",y.ngrid,y.ngrid)
+            ones = Array{Float64,1}(undef,y.n)
+            result = Array{Float64,1}(undef,y.n)
+            
+            ones .= 1.0
+            indefinite_integral!(result, ones, y, y_spectral)
+            #println(result)
+            if y.name=="vpa" || y.name =="z"
+                yllim = y.grid[1]
+            else
+                yllim = 0.0
+            end
+            @. result += yllim
+            #println(result)
+            @. result -= y.grid
+            println("max(abs(err)) integration of constant: ",maximum(abs.(result)))
+
+            cosine = Array{Float64,1}(undef,y.n)
+            for iy in 1:y.n
+                cosine[iy] = cos(y.grid[iy]-yllim)
+            end
+            indefinite_integral!(result, cosine, y, y_spectral)
+            for iy in 1:y.n
+                result[iy] -= sin(y.grid[iy]-yllim)
+            end
+            println("max(abs(err)) integration of cosine: ",maximum(abs.(result)))
+            #println(result)
             #print_matrix(y_spectral.radau.M0,"local radau mass matrix M0",y.ngrid,y.ngrid)
             #print_matrix(y_spectral.radau.M1,"local radau mass matrix M1",y.ngrid,y.ngrid)
             #print_matrix(y_spectral.lobatto.M0,"local mass matrix M0",y.ngrid,y.ngrid)
diff --git a/test_scripts/phi_from_Epar_test.jl b/test_scripts/phi_from_Epar_test.jl
new file mode 100644
index 000000000..25a066e8a
--- /dev/null
+++ b/test_scripts/phi_from_Epar_test.jl
@@ -0,0 +1,134 @@
+export test_phi_from_Epar
+
+using Printf
+using Plots
+using LaTeXStrings
+using MPI
+using Measures
+using Dates
+
+import moment_kinetics
+using moment_kinetics.array_allocation: allocate_float, allocate_shared_float
+using moment_kinetics.coordinates: define_coordinate, get_coordinate_input
+using moment_kinetics.communication: setup_distributed_memory_MPI, block_rank, block_size, finalize_comms!
+using moment_kinetics.type_definitions: OptionsDict
+using moment_kinetics.looping
+using moment_kinetics.input_structs: set_defaults_and_check_section!, em_fields_input
+using moment_kinetics.em_fields: setup_em_fields, calculate_phi_from_Epar!
+
+function test_phi_from_Epar(;ngrid_z=5, nelement_local_z=1, nelement_global_z=1, element_spacing_option_z="uniform", Lz = 1.0, 
+                             ngrid_r=1, nelement_local_r=1, nelement_global_r=1, element_spacing_option_r="uniform", Lr = 1.0,
+                             discretization="chebyshev_pseudospectral", ignore_MPI=false)
+    input_dict = OptionsDict(
+            "r"=>OptionsDict("ngrid"=>ngrid_r, "nelement"=>nelement_global_r,
+                                 "nelement_local"=>nelement_local_r, "L"=>Lr,
+                                 "discretization"=>discretization,
+                                 "element_spacing_option"=>element_spacing_option_r),
+            "z"=>OptionsDict("ngrid"=>ngrid_z, "nelement"=>nelement_global_z,
+                               "nelement_local"=>nelement_local_z, "L"=>Lz,
+                               "discretization"=>discretization,
+                               "element_spacing_option"=>element_spacing_option_z),
+        )
+    
+    # set up distributed-memory MPI information for z and r coords
+    # need grid and MPI information to determine these values 
+    # MRH just put dummy values now 
+    r_coord_input = get_coordinate_input(input_dict, "r"; ignore_MPI=ignore_MPI)
+    z_coord_input = get_coordinate_input(input_dict, "z"; ignore_MPI=ignore_MPI)
+    if ignore_MPI
+        irank_z = irank_r = 0
+        nrank_z = nrank_r = 1
+        comm_sub_z = comm_sub_r = MPI.COMM_NULL
+    else
+        irank_z, nrank_z, comm_sub_z, irank_r, nrank_r, comm_sub_r =
+            setup_distributed_memory_MPI(z_coord_input.nelement,
+                                         z_coord_input.nelement_local,
+                                         r_coord_input.nelement,
+                                         r_coord_input.nelement_local)
+    end
+    z, z_spectral = define_coordinate(z_coord_input; 
+                                      run_directory="", ignore_MPI=ignore_MPI,
+                                      irank=irank_z, nrank=nrank_z, comm=comm_sub_z)
+    r, r_spectral = define_coordinate(r_coord_input;
+                                      run_directory="",ignore_MPI=ignore_MPI,
+                                      irank=irank_r, nrank=nrank_r, comm=comm_sub_r)
+    
+    em_input = set_defaults_and_check_section!(
+                                        input_dict, em_fields_input, "em_fields"
+                                       )
+    looping.setup_loop_ranges!(block_rank[], block_size[];
+                                      s=1,
+                                      sn=1,
+                                      r=r.n, z=z.n, vperp=1, vpa=1,
+                                      vzeta=1, vr=1, vz=1)
+
+    # create the "fields" structure that contains arrays
+    # for the electrostatic potential phi and the electromagnetic fields
+    # set vperp.n = n_ion_species = 1
+    fields = setup_em_fields(1, z.n, r.n, 1,
+                             em_input)
+
+    # setup test
+    phi = fields.phi
+    phi_exact = deepcopy(phi)
+    phi_error = deepcopy(phi)
+    Epar = fields.Ez
+    epsilon = 0.1
+    offset = -3.0
+    begin_serial_region()
+
+    @serial_region begin
+        @loop_r_z ir iz begin
+            zplus = 0.5 + z.grid[iz]/z.L + epsilon 
+            zminus = 0.5 - z.grid[iz]/z.L + epsilon 
+            phi_exact[iz,ir] = sqrt(zplus*zminus) + offset
+            Epar[iz,ir] = -(0.5/z.L)*( sqrt(zminus/zplus) - sqrt(zplus/zminus) ) 
+        end
+    end
+
+    # now calculate phi
+    # calculate_phi_from_Epar!() assumes that phi[1,:] and phi[end,:] have
+    # been set with the correct values by a kinetic boundary condition
+    if z.irank == 0
+        @serial_region begin
+            @loop_r ir begin
+                phi[1,ir] = phi_exact[1,ir]
+        
+            end
+        end
+    end
+    if z.irank == z.nrank - 1
+        @serial_region begin
+            @loop_r ir begin
+                phi[end,ir] = phi_exact[end,ir]
+            end
+        end
+    end
+    calculate_phi_from_Epar!(phi, Epar, r, z, z_spectral)
+
+    # evaluate the error
+    @serial_region begin
+        @loop_r_z ir iz begin
+            phi_error[iz,ir] = phi[iz,ir] - phi_exact[iz,ir]
+        end
+        #println("phi:", phi)
+        #println("phi_exact:", phi_exact)
+        #println("phi_error:", phi_error)
+        println("z.irank: ",z.irank," max(abs(phi_error)): ",maximum(abs.(phi_error) ))
+
+    end
+
+    finalize_comms!()
+end
+# run test by, e.g.,
+# mpirun -n 8 --output-filename testpath julia --project -O3 test_scripts/phi_from_Epar_test.jl
+if abspath(PROGRAM_FILE) == @__FILE__
+    using Pkg
+    Pkg.activate(".")
+    # for precompilation
+    # test_phi_from_Epar()
+    test_phi_from_Epar(ngrid_z=5, nelement_global_z=16, nelement_local_z=2)
+    test_phi_from_Epar(ngrid_z=5, nelement_global_z=16, nelement_local_z=4,
+                       ngrid_r=5, nelement_global_r=2, nelement_local_r=1)
+    
+end
\ No newline at end of file