diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index df4358220..6b127c0b5 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -7,6 +7,8 @@ set(fppFiles
     stdlib_experimental_stats.fypp
     stdlib_experimental_stats_mean.fypp
     stdlib_experimental_stats_var.fypp
+    stdlib_experimental_quadrature.fypp
+    stdlib_experimental_quadrature_trapz.fypp
 )
 
 
diff --git a/src/stdlib_experimental_quadrature.fypp b/src/stdlib_experimental_quadrature.fypp
new file mode 100644
index 000000000..ecb15840f
--- /dev/null
+++ b/src/stdlib_experimental_quadrature.fypp
@@ -0,0 +1,75 @@
+#:include "common.fypp"
+
+module stdlib_experimental_quadrature
+    use stdlib_experimental_kinds, only: sp, dp, qp
+
+    implicit none
+
+    private
+
+    ! array integration
+    public :: trapz
+    public :: trapz_weights
+    public :: simps
+    public :: simps_weights
+
+
+    interface trapz
+        #:for KIND in REAL_KINDS
+        pure module function trapz_dx_${KIND}$(y, dx) result(integral)
+            real(${KIND}$), dimension(:), intent(in) :: y
+            real(${KIND}$), intent(in) :: dx
+            real(${KIND}$) :: integral
+        end function trapz_dx_${KIND}$
+        #:endfor
+        #:for KIND in REAL_KINDS
+        pure module function trapz_x_${KIND}$(y, x) result(integral)
+            real(${KIND}$), dimension(:), intent(in) :: y
+            real(${KIND}$), dimension(size(y)), intent(in) :: x
+            real(${KIND}$) :: integral
+        end function trapz_x_${KIND}$
+        #:endfor
+    end interface trapz
+
+
+    interface trapz_weights
+        #:for KIND in REAL_KINDS
+        pure module function trapz_weights_${KIND}$(x) result(w)
+            real(${KIND}$), dimension(:), intent(in) :: x
+            real(${KIND}$), dimension(size(x)) :: w
+        end function trapz_weights_${KIND}$
+        #:endfor
+    end interface trapz_weights
+
+
+    interface simps
+        #:for KIND in REAL_KINDS
+        pure module function simps_dx_${KIND}$(y, dx, even) result(integral)
+            real(${KIND}$), dimension(:), intent(in) :: y
+            real(${KIND}$), intent(in) :: dx
+            integer, intent(in), optional :: even
+            real(${KIND}$) :: integral
+        end function simps_dx_${KIND}$
+        #:endfor
+        #:for KIND in REAL_KINDS
+        pure module function simps_x_${KIND}$(y, x, even) result(integral)
+            real(${KIND}$), dimension(:), intent(in) :: y
+            real(${KIND}$), dimension(size(y)), intent(in) :: x
+            integer, intent(in), optional :: even
+            real(${KIND}$) :: integral
+        end function simps_x_${KIND}$
+        #:endfor
+    end interface simps
+
+
+    interface simps_weights
+        #:for KIND in REAL_KINDS
+        pure module function simps_weights_${KIND}$(x, even) result(w)
+            real(${KIND}$), dimension(:), intent(in) :: x
+            real(${KIND}$), dimension(size(x)) :: w
+            integer, intent(in), optional :: even
+        end function simps_weights_${KIND}$
+        #:endfor
+    end interface simps_weights
+
+end module stdlib_experimental_quadrature
diff --git a/src/stdlib_experimental_quadrature.md b/src/stdlib_experimental_quadrature.md
new file mode 100644
index 000000000..c62638f64
--- /dev/null
+++ b/src/stdlib_experimental_quadrature.md
@@ -0,0 +1,141 @@
+# Numerical integration
+
+## Implemented
+
+* `trapz`
+* `trapz_weights`
+
+## `trapz` - integrate sampled values using trapezoidal rule
+
+Returns the trapezoidal rule integral of an array `y` representing discrete samples of a function. The integral is computed assuming either equidistant abscissas with spacing `dx` or arbitary abscissas `x`.
+
+### Syntax
+
+`result = trapz(y, x)`
+
+`result = trapz(y, dx)`
+
+### Arguments
+
+`y`: Shall be a rank-one array of type `real`.
+
+`x`: Shall be a rank-one array of type `real` having the same kind and size as `y`. 
+
+`dx`: Shall be a scalar of type `real` having the same kind as `y`.
+
+### Return value
+
+The result is a scalar of type `real` having the same kind as `y`.
+
+If the size of `y` is zero or one, the result is zero.
+
+### Example
+
+```fortran
+program demo_trapz
+    use stdlib_experimental_quadrature, only: trapz
+    implicit none
+    real :: x(5) = [0., 1., 2., 3., 4.]
+    real :: y(5) = x**2
+    print *, trapz(y, x) 
+! 22.0
+    print *, trapz(y, 0.5) 
+! 11.0
+end program
+```
+
+## `trapz_weights` - trapezoidal rule weights for given abscissas
+
+Given an array of abscissas `x`, computes the array of weights `w` such that if `y` represented function values tabulated at `x`, then `sum(w*y)` produces a trapezoidal rule approximation to the integral.
+
+### Syntax
+
+`result = trapz_weights(x)`
+
+### Arguments
+
+`x`: Shall be a rank-one array of type `real`.
+
+### Return value
+
+The result is a `real` array with the same size and kind as `x`.
+
+If the size of `x` is one, then the sole element of the result is zero.
+
+### Example
+
+```fortran
+program demo_trapz_weights
+    use stdlib_experimental_quadrature, only: trapz_weights
+    implicit none
+    real :: x(5) = [0., 1., 2., 3., 4.]
+    real :: y(5) = x**2
+    real :: w(5) 
+    w = trapz_weight(x)
+    print *, dot_product(w, y)
+! 22.0
+end program
+
+```
+
+# `simps` - integrate sampled values using Simpson's rule
+
+Returns the Simpson's rule integral of an array `y` representing discrete samples of a function. The integral is computed assuming either equidistant abscissas with spacing `dx` or arbitary abscissas `x`. 
+
+Simpson's rule is defined for odd-length arrays only. For even-length arrays, an optional argument `even` may be used to specify at which index to replace Simpson's rule with Simpson's 3/8 rule. The 3/8 rule will be used for the array section `y(even:even+4)` and the ordinary Simpon's rule will be used elsewhere.
+
+### Syntax
+
+`result = simps(y, x [, even])`
+
+`result = simps(y, dx [, even])`
+
+### Arguments
+
+`y`: Shall be a rank-one array of type `real`.
+
+`x`: Shall be a rank-one array of type `real` having the same kind and size as `y`. 
+
+`dx`: Shall be a scalar of type `real` having the same kind as `y`.
+
+`even`: (Optional) Shall be a scalar integer of default kind. Its default value is `1`.
+
+### Return value
+
+The result is a scalar of type `real` having the same kind as `y`.
+
+If the size of `y` is zero or one, the result is zero.
+
+If the size of `y` is two, the result is the same as if `trapz` had been called instead, regardless of the value of `even`.
+
+### Example
+
+TBD
+
+# `simps_weights` - Simpson's rule weights for given abscissas
+
+Given an array of abscissas `x`, computes the array of weights `w` such that if `y` represented function values tabulated at `x`, then `sum(w*y)` produces a Simpson's rule approximation to the integral.
+
+Simpson's rule is defined for odd-length arrays only. For even-length arrays, an optional argument `even` may be used to specify at which index to replace Simpson's rule with Simpson's 3/8 rule. The 3/8 rule will be used for the array section `x(even:even+4)` and the ordinary Simpon's rule will be used elsewhere.
+
+### Syntax
+
+`result = simps_weights(x [, even])`
+
+### Arguments
+
+`x`: Shall be a rank-one array of type `real`.
+
+`even`: (Optional) Shall be a scalar integer of default kind. Its default value is `1`.
+
+### Return value
+
+The result is a `real` array with the same size and kind as `x`.
+
+If the size of `x` is one, then the sole element of the result is zero.
+
+If the size of `x` is two, then the result is the same as if `trapz_weights` had been called instead.
+
+### Example
+
+TBD
diff --git a/src/stdlib_experimental_quadrature_trapz.fypp b/src/stdlib_experimental_quadrature_trapz.fypp
new file mode 100644
index 000000000..0ea66fe5e
--- /dev/null
+++ b/src/stdlib_experimental_quadrature_trapz.fypp
@@ -0,0 +1,86 @@
+#:include "common.fypp"
+
+submodule (stdlib_experimental_quadrature) stdlib_experimental_quadrature_trapz
+
+    implicit none
+
+contains
+
+    #:for KIND in REAL_KINDS
+    
+    pure module function trapz_dx_${KIND}$(y, dx) result(integral)
+        real(${KIND}$), dimension(:), intent(in) :: y
+        real(${KIND}$), intent(in) :: dx
+        real(${KIND}$) :: integral
+
+        integer :: n
+
+        n = size(y)
+        
+        select case (n)
+        case (0:1)
+            integral = 0.0_${KIND}$
+        case (2)
+            integral = 0.5_${KIND}$*dx*(y(1) + y(2))
+        case default
+            integral = dx*(sum(y(2:n-1)) + 0.5_${KIND}$*(y(1) + y(n)))
+        end select
+    end function trapz_dx_${KIND}$
+    
+    #:endfor
+    #:for KIND in REAL_KINDS
+    
+    pure module function trapz_x_${KIND}$(y, x) result(integral)
+        real(${KIND}$), dimension(:), intent(in) :: y
+        real(${KIND}$), dimension(size(y)), intent(in) :: x
+        real(${KIND}$) :: integral
+
+        integer :: i
+        integer :: n
+
+        n = size(y)
+
+        select case (n)
+        case (0:1)
+            integral = 0.0_${KIND}$
+        case (2)
+            integral = 0.5_${KIND}$*(x(2) - x(1))*(y(1) + y(2))
+        case default
+            integral = 0.0_${KIND}$
+            do i = 2, n
+                integral = integral + (x(i) - x(i-1))*(y(i) + y(i-1))
+            end do
+            integral = 0.5_${KIND}$*integral
+        end select
+    end function trapz_x_${KIND}$
+
+    #:endfor
+    #:for KIND in REAL_KINDS
+
+    pure module function trapz_weights_${KIND}$(x) result(w)
+        real(${KIND}$), dimension(:), intent(in) :: x
+        real(${KIND}$), dimension(size(x)) :: w
+
+        integer :: i
+        integer :: n
+
+        n = size(x)
+
+        select case (n)
+        case (0)
+            ! no action needed
+        case (1)
+            w(1) = 0.0_${KIND}$
+        case (2)
+            w = 0.5_${KIND}$*(x(2) - x(1))
+        case default
+            w(1) = 0.5_${KIND}$*(x(2) - x(1))
+            w(n) = 0.5_${KIND}$*(x(n) - x(n-1))
+            do i = 2, size(x)-1
+                w(i) = 0.5_${KIND}$*(x(i+1) - x(i-1))
+            end do
+        end select
+    end function trapz_weights_${KIND}$
+
+#:endfor
+end submodule stdlib_experimental_quadrature_trapz
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
index df5bd0a09..f3b7d434c 100644
--- a/src/tests/CMakeLists.txt
+++ b/src/tests/CMakeLists.txt
@@ -11,6 +11,7 @@ add_subdirectory(io)
 add_subdirectory(optval)
 add_subdirectory(stats)
 add_subdirectory(system)
+add_subdirectory(quadrature)
 
 ADDTEST(always_skip)
 set_tests_properties(always_skip PROPERTIES SKIP_RETURN_CODE 77)
diff --git a/src/tests/quadrature/CMakeLists.txt b/src/tests/quadrature/CMakeLists.txt
new file mode 100644
index 000000000..08f1ce4c6
--- /dev/null
+++ b/src/tests/quadrature/CMakeLists.txt
@@ -0,0 +1,2 @@
+ADDTEST(trapz)
+
diff --git a/src/tests/quadrature/Makefile.manual b/src/tests/quadrature/Makefile.manual
new file mode 100644
index 000000000..83061a7bf
--- /dev/null
+++ b/src/tests/quadrature/Makefile.manual
@@ -0,0 +1,3 @@
+PROGS_SRC = test_trapz.f90
+
+include ../Makefile.manual.test.mk
diff --git a/src/tests/quadrature/test_trapz.f90 b/src/tests/quadrature/test_trapz.f90
new file mode 100644
index 000000000..1265e4b6b
--- /dev/null
+++ b/src/tests/quadrature/test_trapz.f90
@@ -0,0 +1,236 @@
+program test_trapz
+    use stdlib_experimental_kinds, only: sp, dp, qp
+    use stdlib_experimental_error, only: assert
+    use stdlib_experimental_quadrature, only: trapz, trapz_weights
+
+    implicit none
+
+    call test_trapz_sp
+    call test_trapz_dp
+    call test_trapz_qp
+
+    call test_trapz_weights_sp
+    call test_trapz_weights_dp
+    call test_trapz_weights_qp
+
+    call test_trapz_zero_sp
+    call test_trapz_zero_dp
+    call test_trapz_zero_qp
+
+contains
+
+    subroutine test_trapz_sp
+        integer, parameter :: n = 17
+        real(sp), dimension(n) :: y
+        real(sp), dimension(n) :: x
+        real(sp) :: val
+        real(sp) :: ans
+        integer :: i
+
+        print *, "test_trapz_sp"
+
+        y = [(real(i-1, sp), i = 1, n)]
+
+        val = trapz(y, 1.0_sp)
+        ans = 128.0_sp
+        call assert(abs(val - ans) < epsilon(ans))
+
+        val = trapz(y, 0.5_sp)
+        ans = 64.0_sp
+        call assert(abs(val - ans) < epsilon(ans))
+
+        x = [((i-1)*4.0_sp/real(n-1, sp), i = 1, n)]
+        val = trapz(y, x)
+        ans = 32.0_sp
+        call assert(abs(val - ans) < epsilon(ans))
+
+        x = y**2
+        val = trapz(y, x)
+        ans = 2728.0_sp
+        call assert(abs(val - ans) < epsilon(ans))
+    end subroutine test_trapz_sp
+
+    subroutine test_trapz_dp
+        integer, parameter :: n = 17
+        real(dp), dimension(n) :: y
+        real(dp), dimension(n) :: x
+        real(dp) :: val
+        real(dp) :: ans
+        integer :: i
+
+        print *, "test_trapz_dp"
+
+        y = [(real(i-1, dp), i = 1, n)]
+
+        val = trapz(y, 1.0_dp)
+        ans = 128.0_dp
+        call assert(abs(val - ans) < epsilon(ans))
+
+        val = trapz(y, 0.5_dp)
+        ans = 64.0_dp
+        call assert(abs(val - ans) < epsilon(ans))
+
+        x = [((i-1)*4.0_dp/real(n-1, dp), i = 1, n)]
+        val = trapz(y, x)
+        ans = 32.0_dp
+        call assert(abs(val - ans) < epsilon(ans))
+
+        x = y**2
+        val = trapz(y, x)
+        ans = 2728.0_sp
+        call assert(abs(val - ans) < epsilon(ans))
+    end subroutine test_trapz_dp
+
+
+    subroutine test_trapz_qp
+        integer, parameter :: n = 17
+        real(qp), dimension(n) :: y
+        real(qp), dimension(n) :: x
+        real(qp) :: val
+        real(qp) :: ans
+        integer :: i
+
+        print *, "test_trapz_qp"
+
+        y = [(real(i-1, qp), i = 1, n)]
+
+        val = trapz(y, 1.0_qp)
+        ans = 128.0_qp
+        call assert(abs(val - ans) < epsilon(ans))
+
+        val = trapz(y, 0.5_qp)
+        ans = 64.0_qp
+        call assert(abs(val - ans) < epsilon(ans))
+
+        x = [((i-1)*4.0_qp/real(n-1, qp), i = 1, n)]
+        val = trapz(y, x)
+        ans = 32.0_qp
+        call assert(abs(val - ans) < epsilon(ans))
+
+        x = y**2
+        val = trapz(y, x)
+        ans = 2728.0_qp
+        call assert(abs(val - ans) < epsilon(ans))
+    end subroutine test_trapz_qp
+
+
+    subroutine test_trapz_weights_sp
+        integer, parameter :: n = 17
+        real(sp), dimension(n) :: y
+        real(sp), dimension(n) :: x
+        real(sp), dimension(n) :: w
+        integer :: i
+        real(sp) :: val
+        real(sp) :: ans
+
+        print *, "test_trapz_weights_sp"
+
+        y = [(real(i-1, sp), i = 1, n)]
+
+        x = y
+        w = trapz_weights(x)
+        val = dot_product(w, y)
+        ans = trapz(y, x)
+        call assert(abs(val - ans) < epsilon(ans))
+
+        x = y**2
+        w = trapz_weights(x)
+        val = dot_product(w, y)
+        ans = trapz(y, x)
+        call assert(abs(val - ans) < epsilon(ans))
+
+    end subroutine test_trapz_weights_sp
+
+
+    subroutine test_trapz_weights_dp
+        integer, parameter :: n = 17
+        real(dp), dimension(n) :: y
+        real(dp), dimension(n) :: x
+        real(dp), dimension(n) :: w
+        integer :: i
+        real(dp) :: val
+        real(dp) :: ans
+
+        print *, "test_trapz_weights_dp"
+
+        y = [(real(i-1, dp), i = 1, n)]
+
+        x = y
+        w = trapz_weights(x)
+        val = dot_product(w, y)
+        ans = trapz(y, x)
+        call assert(abs(val - ans) < epsilon(ans))
+
+        x = y**2
+        w = trapz_weights(x)
+        val = dot_product(w, y)
+        ans = trapz(y, x)
+        call assert(abs(val - ans) < epsilon(ans))
+
+    end subroutine test_trapz_weights_dp
+
+
+    subroutine test_trapz_weights_qp
+        integer, parameter :: n = 17
+        real(qp), dimension(n) :: y
+        real(qp), dimension(n) :: x
+        real(qp), dimension(n) :: w
+        integer :: i
+        real(qp) :: val
+        real(qp) :: ans
+
+        print *, "test_trapz_weights_qp"
+
+        y = [(real(i-1, qp), i = 1, n)]
+
+        x = y
+        w = trapz_weights(x)
+        val = dot_product(w, y)
+        ans = trapz(y, x)
+        call assert(abs(val - ans) < epsilon(ans))
+
+        x = y**2
+        w = trapz_weights(x)
+        val = dot_product(w, y)
+        ans = trapz(y, x)
+        call assert(abs(val - ans) < epsilon(ans))
+
+    end subroutine test_trapz_weights_qp
+
+
+    subroutine test_trapz_zero_sp
+        real(sp), dimension(0) :: a
+
+        print *, "test_trapz_zero_sp"
+
+        call assert(abs(trapz(a, 1.0_sp)) < epsilon(0.0_sp))
+        call assert(abs(trapz([1.0_sp], 1.0_sp)) < epsilon(0.0_sp))
+        call assert(abs(trapz(a, a)) < epsilon(0.0_sp))
+        call assert(abs(trapz([1.0_sp], [1.0_sp])) < epsilon(0.0_sp))
+    end subroutine test_trapz_zero_sp
+
+
+    subroutine test_trapz_zero_dp
+        real(dp), dimension(0) :: a
+
+        print *, "test_trapz_zero_dp"
+
+        call assert(abs(trapz(a, 1.0_dp)) < epsilon(0.0_dp))
+        call assert(abs(trapz([1.0_dp], 1.0_dp)) < epsilon(0.0_dp))
+        call assert(abs(trapz(a, a)) < epsilon(0.0_dp))
+        call assert(abs(trapz([1.0_dp], [1.0_dp])) < epsilon(0.0_dp))
+    end subroutine test_trapz_zero_dp
+
+
+    subroutine test_trapz_zero_qp
+        real(qp), dimension(0) :: a
+
+        print *, "test_trapz_zero_qp"
+
+        call assert(abs(trapz(a, 1.0_qp)) < epsilon(0.0_qp))
+        call assert(abs(trapz([1.0_qp], 1.0_qp)) < epsilon(0.0_qp))
+        call assert(abs(trapz(a, a)) < epsilon(0.0_qp))
+        call assert(abs(trapz([1.0_qp], [1.0_qp])) < epsilon(0.0_qp))
+    end subroutine test_trapz_zero_qp
+
+end program test_trapz