diff --git a/doc/specs/stdlib_experimental_quadrature.md b/doc/specs/stdlib_experimental_quadrature.md
index 92a35849f..9bde409ce 100644
--- a/doc/specs/stdlib_experimental_quadrature.md
+++ b/doc/specs/stdlib_experimental_quadrature.md
@@ -76,26 +76,26 @@ program demo_trapz_weights
     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)
+    w = trapz_weights(x)
+    print *, sum(w*y)
 ! 22.0
 end program demo_trapz_weights
 
 ```
 
-## `simps` - integrate sampled values using Simpson's rule (to be implemented)
+## `simps` - integrate sampled values using Simpson's rule
 
 ### Description
 
 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.
+Simpson's ordinary ("1/3") rule is used for odd-length arrays. For even-length arrays, Simpson's 3/8 rule is also utilized in a way that depends on the value of `even`. If `even` is negative (positive), the 3/8 rule is used at the beginning (end) of the array. If `even` is zero or not present, the result is as if the 3/8 rule were first used at the beginning of the array, then at the end of the array, and these two results were averaged.
 
 ### Syntax
 
-`result = simps(y, x [, even])`
+`result = [[stdlib_experimental_quadrature(module):simps(interface)]](y, x [, even])`
 
-`result = simps(y, dx [, even])`
+`result = [[stdlib_experimental_quadrature(module):simps(interface)]](y, dx [, even])`
 
 ### Arguments
 
@@ -105,7 +105,7 @@ Simpson's rule is defined for odd-length arrays only. For even-length arrays, an
 
 `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`.
+`even`: (Optional) Shall be a default-kind `integer`.
 
 ### Return value
 
@@ -113,29 +113,40 @@ 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`.
+If the size of `y` is two, the result is the same as if `trapz` had been called instead.
 
 ### Example
 
-TBD
+```fortran
+program demo_simps
+    use stdlib_experimental_quadrature, only: simps
+    implicit none
+    real :: x(5) = [0., 1., 2., 3., 4.]
+    real :: y(5) = 3.*x**2
+    print *, simps(y, x) 
+! 64.0
+    print *, simps(y, 0.5) 
+! 32.0
+end program demo_simps
+```
 
-## `simps_weights` - Simpson's rule weights for given abscissas (to be implemented)
+## `simps_weights` - Simpson's rule weights for given abscissas
 
 ### Description
 
 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.
+Simpson's ordinary ("1/3") rule is used for odd-length arrays. For even-length arrays, Simpson's 3/8 rule is also utilized in a way that depends on the value of `even`. If `even` is negative (positive), the 3/8 rule is used at the beginning (end) of the array and the 1/3 rule used elsewhere. If `even` is zero or not present, the result is as if the 3/8 rule were first used at the beginning of the array, then at the end of the array, and then these two results were averaged.
 
 ### Syntax
 
-`result = simps_weights(x [, even])`
+`result = [[stdlib_experimental_quadrature(module):simps_weights(interface)]](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`.
+`even`: (Optional) Shall be a default-kind `integer`.
 
 ### Return value
 
@@ -147,4 +158,15 @@ If the size of `x` is two, then the result is the same as if `trapz_weights` had
 
 ### Example
 
-TBD
+```fortran
+program demo_simps_weights
+    use stdlib_experimental_quadrature, only: simps_weights
+    implicit none
+    real :: x(5) = [0., 1., 2., 3., 4.]
+    real :: y(5) = 3.*x**2
+    real :: w(5) 
+    w = simps_weights(x)
+    print *, sum(w*y)
+! 64.0
+end program demo_simps_weights
+```
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 51fd01d62..c4cff0cfd 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -13,6 +13,7 @@ set(fppFiles
     stdlib_experimental_stats_var.fypp
     stdlib_experimental_quadrature.fypp
     stdlib_experimental_quadrature_trapz.fypp
+    stdlib_experimental_quadrature_simps.fypp
 )
 
 
diff --git a/src/stdlib_experimental_quadrature.fypp b/src/stdlib_experimental_quadrature.fypp
index 075c19b11..42270dfe8 100644
--- a/src/stdlib_experimental_quadrature.fypp
+++ b/src/stdlib_experimental_quadrature.fypp
@@ -1,5 +1,4 @@
 #:include "common.fypp"
-
 module stdlib_experimental_quadrature
     !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description))
     use stdlib_experimental_kinds, only: sp, dp, qp
@@ -18,19 +17,19 @@ module stdlib_experimental_quadrature
     interface trapz
         !! Integrates sampled values using trapezoidal rule
         !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description))
-        #: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}$
+        #:for k1, t1 in REAL_KINDS_TYPES
+        pure module function trapz_dx_${k1}$(y, dx) result(integral)
+          ${t1}$, dimension(:), intent(in) :: y
+          ${t1}$, intent(in) :: dx
+          ${t1}$ :: integral
+        end function trapz_dx_${k1}$      
         #: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}$
+        #:for k1, t1 in REAL_KINDS_TYPES
+        module function trapz_x_${k1}$(y, x) result(integral)
+            ${t1}$, dimension(:), intent(in) :: y
+            ${t1}$, dimension(:), intent(in) :: x
+            ${t1}$ :: integral
+        end function trapz_x_${k1}$
         #:endfor
     end interface trapz
 
@@ -38,43 +37,62 @@ module stdlib_experimental_quadrature
     interface trapz_weights
         !! Integrates sampled values using trapezoidal rule weights for given abscissas
         !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description_1))
-        #: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}$
+        #:for k1, t1 in REAL_KINDS_TYPES
+        pure module function trapz_weights_${k1}$(x) result(w)
+            ${t1}$, dimension(:), intent(in) :: x
+            ${t1}$, dimension(size(x)) :: w
+        end function trapz_weights_${k1}$
         #: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
+        !! Integrates sampled values using Simpson's rule
+        !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description_3))
+        ! "recursive" is an implementation detail
+        #:for k1, t1 in REAL_KINDS_TYPES
+        pure recursive module function simps_dx_${k1}$(y, dx, even) result(integral)
+            ${t1}$, dimension(:), intent(in) :: y
+            ${t1}$, intent(in) :: dx
             integer, intent(in), optional :: even
-            real(${KIND}$) :: integral
-        end function simps_dx_${KIND}$
+            ${t1}$ :: integral
+        end function simps_dx_${k1}$
         #: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
+        #:for k1, t1 in REAL_KINDS_TYPES
+        recursive module function simps_x_${k1}$(y, x, even) result(integral)
+            ${t1}$, dimension(:), intent(in) :: y
+            ${t1}$, dimension(:), intent(in) :: x
             integer, intent(in), optional :: even
-            real(${KIND}$) :: integral
-        end function simps_x_${KIND}$
+            ${t1}$ :: integral
+        end function simps_x_${k1}$
         #: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
+        !! Integrates sampled values using trapezoidal rule weights for given abscissas
+        !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description_3))
+        #:for k1, t1 in REAL_KINDS_TYPES
+        pure recursive module function simps_weights_${k1}$(x, even) result(w)
+            ${t1}$, dimension(:), intent(in) :: x
             integer, intent(in), optional :: even
-        end function simps_weights_${KIND}$
+            ${t1}$, dimension(size(x)) :: w
+        end function simps_weights_${k1}$
         #:endfor
     end interface simps_weights
 
+
+    ! Interface for a simple f(x)-style integrand function.
+    ! Could become fancier as we learn about the performance
+    ! ramifications of different ways to do callbacks.
+    abstract interface
+        #:for k1, t1 in REAL_KINDS_TYPES
+        pure function integrand_${k1}$(x) result(f)
+            import :: ${k1}$
+            ${t1}$, intent(in) :: x
+            ${t1}$ :: f
+        end function integrand_${k1}$
+        #:endfor
+    end interface
+
 end module stdlib_experimental_quadrature
diff --git a/src/stdlib_experimental_quadrature_simps.fypp b/src/stdlib_experimental_quadrature_simps.fypp
new file mode 100644
index 000000000..8ff4c5437
--- /dev/null
+++ b/src/stdlib_experimental_quadrature_simps.fypp
@@ -0,0 +1,278 @@
+#:include "common.fypp"
+
+submodule (stdlib_experimental_quadrature) stdlib_experimental_quadrature_simps
+    use stdlib_experimental_error, only: check
+    implicit none
+
+    ! internal use only
+    interface simps38
+    #:for k1, t1 in REAL_KINDS_TYPES
+        module procedure simps38_dx_${k1}$
+        module procedure simps38_x_${k1}$
+    #:endfor
+    end interface simps38
+
+    ! internal use only
+    interface simps38_weights
+    #:for k1, t1 in REAL_KINDS_TYPES
+        module procedure simps38_weights_${k1}$
+    #:endfor
+    end interface simps38_weights
+
+contains
+
+#:for k1, t1 in REAL_KINDS_TYPES
+
+    pure recursive module function simps_dx_${k1}$(y, dx, even) result(integral)
+        ${t1}$, dimension(:), intent(in) :: y
+        ${t1}$, intent(in) :: dx
+        integer, intent(in), optional :: even
+        ${t1}$ :: integral
+
+        integer :: n
+
+        n = size(y)
+
+        select case (n)
+        case (0:1)
+            integral = 0.0_${k1}$
+        case (2)
+            integral = 0.5_${k1}$*dx*(y(1) + y(2))
+        case (3)
+            integral = dx/3.0_${k1}$*(y(1) + 4*y(2) + y(3))
+        case (4)
+            integral = simps38(y, dx)
+        ! case (5) not needed; handled by default
+        case (6) ! needs special handling because of averaged 3/8's rule case
+            if (present(even)) then
+                if (even < 0) then
+                    ! 3/8 rule on left
+                    integral =  simps38(y(1:4), dx) + simps(y(4:6), dx)
+                    return
+                else if (even > 0) then
+                    ! 3/8 rule on right
+                    integral = simps(y(1:3), dx) + simps38(y(3:6), dx)
+                    return
+                else
+                    ! fall through
+                end if
+            end if
+            ! either `even` not present or is zero
+            ! equivalent to averaging left and right
+            integral = dx/48.0_${k1}$ * (17*(y(1) + y(6)) + 59*(y(2) + y(5)) + 44*(y(3) + y(4)))
+        case default
+            if (mod(n, 2) == 1) then
+                integral = dx/3.0_${k1}$*(y(1) + 4*sum(y(2:n-1:2)) + 2*sum(y(3:n-2:2)) + y(n))
+            else
+                if (present(even)) then
+                    if (even < 0) then
+                        ! 3/8th rule on left
+                        integral = simps38(y(1:4), dx) + simps(y(4:n), dx)
+                        return
+                    else if (even > 0) then
+                        ! 3/8 rule on right
+                        integral = simps(y(1:n-3), dx) + simps38(y(n-3:n), dx)
+                        return
+                    else
+                        ! fall through
+                    end if
+                end if
+                ! either `even` not present or is zero
+                ! equivalent to averaging left and right
+                integral = dx/48.0_${k1}$ * (17*(y(1) + y(n)) + 59*(y(2) + y(n-1)) &
+                        + 43*(y(3) + y(n-2)) + 49*(y(4) + y(n-3)) + 48*sum(y(5:n-4)))
+            end if
+        end select
+    end function simps_dx_${k1}$
+
+#:endfor
+#:for k1, t1 in REAL_KINDS_TYPES
+
+    recursive module function simps_x_${k1}$(y, x, even) result(integral)
+        ${t1}$, dimension(:), intent(in) :: y
+        ${t1}$, dimension(:), intent(in) :: x
+        integer, intent(in), optional :: even
+        ${t1}$ :: integral
+
+        integer :: i
+        integer :: n
+
+        ${t1}$ :: h1, h2
+        ${t1}$ :: a, b, c
+
+        n = size(y)
+        call check(size(x) == n, "simps: Arguments `x` and `y` must be the same size.")
+
+        select case (n)
+        case (0:1)
+            integral = 0.0_${k1}$
+        case (2)
+            integral = 0.5_${k1}$*(x(2) - x(1))*(y(1) + y(2))
+        case (3)
+            h1 = x(2) - x(1)
+            h2 = x(3) - x(2)
+            a = (2*h1**2 + h1*h2 - h2**2)/(6*h1)
+            b = (h1+h2)**3/(6*h1*h2)
+            c = (2*h2**2 + h1*h2 - h1**2)/(6*h2)
+            integral = a*y(1) + b*y(2) + c*y(3)
+        case (4)
+            integral = simps38(y, x)
+        ! case (6) unneeded; handled by default
+        case default
+            if (mod(n, 2) == 1) then
+                integral = 0.0_${k1}$
+                do i = 1, n-2, 2
+                    h1 = x(i+1) - x(i)
+                    h2 = x(i+2) - x(i+1)
+                    a = (2*h1**2 + h1*h2 - h2**2)/(6*h1)
+                    b = (h1+h2)**3/(6*h1*h2)
+                    c = (2*h2**2 + h1*h2 - h1**2)/(6*h2)
+                    integral = integral + a*y(i) + b*y(i+1) + c*y(i+2)
+                end do
+            else
+                if (present(even)) then
+                    if (even < 0) then
+                        ! 3/8 rule on left
+                        integral = simps38(y(1:4), x(1:4)) + simps(y(4:n), x(4:n))
+                        return
+                    else if (even > 0) then
+                        ! 3/8 rule on right
+                        integral = simps(y(1:n-3), x(1:n-3)) + simps38(y(n-3:n), x(n-3:n))
+                        return
+                    else
+                        ! fall through
+                    end if
+                end if
+                ! either `even` not present or is zero
+                integral = 0.5_${k1}$ * ( simps38(y(1:4), x(1:4)) + simps(y(4:n), x(4:n)) &
+                    + simps(y(1:n-3), x(1:n-3)) + simps38(y(n-3:n), x(n-3:n)) )
+            end if
+        end select
+    end function simps_x_${k1}$
+
+#:endfor
+#:for k1, t1 in REAL_KINDS_TYPES
+
+    pure recursive module function simps_weights_${k1}$(x, even) result(w)
+        ${t1}$, dimension(:), intent(in) :: x
+        integer, intent(in), optional :: even
+        ${t1}$, dimension(size(x)) :: w
+
+        integer :: i, n
+        ${t1}$ :: h1, h2
+
+        n = size(x)
+
+        select case (n)
+        case (0)
+            ! no action needed
+        case (1)
+            w(1) = 0.0_${k1}$
+        case (2)
+            w = 0.5_${k1}$*(x(2) - x(1))
+        case (3)
+            h1 = x(2) - x(1)
+            h2 = x(3) - x(2)
+            w(1) = (2*h1**2 + h1*h2 - h2**2)/(6*h1)
+            w(2) = (h1+h2)**3/(6*h1*h2)
+            w(3) = (2*h2**2 + h1*h2 - h1**2)/(6*h2)
+        case (4)
+            w = simps38_weights(x)
+        case default
+            if (mod(n, 2) == 1) then
+                w = 0.0_${k1}$
+                do i = 1, n-2, 2
+                    h1 = x(i+1) - x(i)
+                    h2 = x(i+2) - x(i+1)
+                    w(i) = w(i) + (2*h1**2 + h1*h2 - h2**2)/(6*h1)
+                    w(i+1) = w(i+1) + (h1+h2)**3/(6*h1*h2)
+                    w(i+2) = w(i+2) + (2*h2**2 + h1*h2 - h1**2)/(6*h2)
+                end do
+            else
+                if (present(even)) then
+                    if (even < 0) then
+                        ! 3/8 rule on left
+                        w = 0.0_${k1}$
+                        w(1:4) = simps38_weights(x(1:4))
+                        w(4:n) = w(4:n) + simps_weights(x(4:n)) ! position 4 needs both rules
+                        return
+                    else if (even > 0) then
+                        ! 3/8 rule on right
+                        w = 0.0_${k1}$
+                        w(1:n-3) = simps_weights(x(1:n-3))
+                        w(n-3:n) = w(n-3:n) + simps38_weights(x(n-3:n)) ! position n-3 needs both rules
+                        return
+                    else
+                        ! fall through
+                    end if
+                end if
+                ! either `even` not present or is zero
+                w = 0.0_${k1}$
+                ! 3/8 rule on left
+                w(1:4) = simps38_weights(x(1:4))
+                w(4:n) = w(4:n) + simps_weights(x(4:n))
+                ! 3/8 rule on right
+                w(1:n-3) = w(1:n-3) + simps_weights(x(1:n-3))
+                w(n-3:n) = w(n-3:n) + simps38_weights(x(n-3:n))
+                ! average
+                w = 0.5_${k1}$ * w
+            end if
+        end select
+    end function simps_weights_${k1}$
+
+#:endfor
+#:for k1, t1 in REAL_KINDS_TYPES
+
+    pure function simps38_dx_${k1}$(y, dx) result (integral)
+        ${t1}$, dimension(4), intent(in) :: y
+        ${t1}$, intent(in) :: dx
+        ${t1}$ :: integral
+
+        integral = 3.0_${k1}$*dx/8.0_${k1}$ * (y(1) + y(4) + 3*(y(2) + y(3)))
+    end function simps38_dx_${k1}$
+
+#:endfor
+#: for k1, t1 in REAL_KINDS_TYPES
+
+    pure function simps38_x_${k1}$(y, x) result(integral)
+        ${t1}$, dimension(4), intent(in) :: y
+        ${t1}$, dimension(4), intent(in) :: x
+        ${t1}$ :: integral
+
+        ${t1}$ :: h1, h2, h3
+        ${t1}$ :: a, b, c, d
+
+        h1 = x(2) - x(1)
+        h2 = x(3) - x(2)
+        h3 = x(4) - x(3)
+
+        a = (h1+h2+h3)*(3*h1**2 + 2*h1*h2 - 2*h1*h3 - h2**2 + h3**2)/(12*h1*(h1+h2))
+        b = (h1+h2-h3)*(h1+h2+h3)**3/(12*h1*h2*(h2+h3))
+        c = (h2+h3-h1)*(h1+h2+h3)**3/(12*h2*h3*(h1+h2))
+        d = (h1+h2+h3)*(3*h3**2 + 2*h2*h3 - 2*h1*h3 - h2**2 + h1**2)/(12*h3*(h2+h3))
+
+        integral = a*y(1) + b*y(2) + c*y(3) + d*y(4)
+    end function simps38_x_${k1}$
+
+#:endfor
+#:for k1, t1 in REAL_KINDS_TYPES
+
+    pure function simps38_weights_${k1}$(x) result(w)
+        ${t1}$, intent(in) :: x(4)
+        ${t1}$ :: w(size(x))
+
+        ${t1}$ :: h1, h2, h3
+
+        h1 = x(2) - x(1)
+        h2 = x(3) - x(2)
+        h3 = x(4) - x(3)
+
+        w(1) = (h1+h2+h3)*(3*h1**2 + 2*h1*h2 - 2*h1*h3 - h2**2 + h3**2)/(12*h1*(h1+h2))
+        w(2) = (h1+h2-h3)*(h1+h2+h3)**3/(12*h1*h2*(h2+h3))
+        w(3) = (h2+h3-h1)*(h1+h2+h3)**3/(12*h2*h3*(h1+h2))
+        w(4) = (h1+h2+h3)*(3*h3**2 + 2*h2*h3 - 2*h1*h3 - h2**2 + h1**2)/(12*h3*(h2+h3))
+    end function simps38_weights_${k1}$
+
+#:endfor
+
+end submodule stdlib_experimental_quadrature_simps
diff --git a/src/stdlib_experimental_quadrature_trapz.fypp b/src/stdlib_experimental_quadrature_trapz.fypp
index 0ea66fe5e..af445b1cb 100644
--- a/src/stdlib_experimental_quadrature_trapz.fypp
+++ b/src/stdlib_experimental_quadrature_trapz.fypp
@@ -1,7 +1,7 @@
 #:include "common.fypp"
 
 submodule (stdlib_experimental_quadrature) stdlib_experimental_quadrature_trapz
-
+    use stdlib_experimental_error, only: check
     implicit none
 
 contains
@@ -29,16 +29,17 @@ contains
     
     #:endfor
     #:for KIND in REAL_KINDS
-    
-    pure module function trapz_x_${KIND}$(y, x) result(integral)
+
+    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}$), dimension(:), intent(in) :: x
         real(${KIND}$) :: integral
 
         integer :: i
         integer :: n
 
         n = size(y)
+        call check(size(x) == n, "trapz: Arguments `x` and `y` must be the same size.")
 
         select case (n)
         case (0:1)
diff --git a/src/tests/quadrature/CMakeLists.txt b/src/tests/quadrature/CMakeLists.txt
index 08f1ce4c6..d8069fb86 100644
--- a/src/tests/quadrature/CMakeLists.txt
+++ b/src/tests/quadrature/CMakeLists.txt
@@ -1,2 +1,2 @@
 ADDTEST(trapz)
-
+ADDTEST(simps)
diff --git a/src/tests/quadrature/test_simps.f90 b/src/tests/quadrature/test_simps.f90
new file mode 100644
index 000000000..b4cdc4287
--- /dev/null
+++ b/src/tests/quadrature/test_simps.f90
@@ -0,0 +1,505 @@
+program test_simps
+    use stdlib_experimental_kinds, only: sp, dp, qp
+    use stdlib_experimental_error, only: check
+    use stdlib_experimental_quadrature, only: simps, simps_weights
+
+    implicit none
+
+    real(sp), parameter :: tol_sp = 1000 * epsilon(1.0_sp)
+    real(dp), parameter :: tol_dp = 1000 * epsilon(1.0_dp)
+    real(qp), parameter :: tol_qp = 1000 *epsilon(1.0_sp)
+
+    call test_simps_sp
+    call test_simps_dp
+    call test_simps_qp
+
+    call test_simps_weights_sp
+    call test_simps_weights_dp
+    call test_simps_weights_qp
+
+    call test_simps_zero_sp
+    call test_simps_zero_dp
+    call test_simps_zero_qp
+
+    call test_simps_even_sp
+    call test_simps_even_dp
+    call test_simps_even_qp
+
+    call test_simps_weights_even_sp
+    call test_simps_weights_even_dp
+    call test_simps_weights_even_qp
+
+    call test_simps_six_sp
+    call test_simps_six_dp
+    call test_simps_six_qp
+
+contains
+
+    subroutine test_simps_sp
+        integer, parameter :: n = 13
+        real(sp), dimension(n) :: y
+        real(sp), dimension(n) :: x
+        real(sp) :: val
+        real(sp) :: ans
+        integer :: i
+
+        print *, "test_simps_sp"
+
+        y = [(real(i-1, sp)**2, i = 1, n)]
+
+        val = simps(y, 1.0_sp)
+        ans = 576.0_sp
+        print *, "  dx=1", val, ans
+        call check(abs(val - ans) < tol_sp)
+
+        val = simps(y, 0.5_sp)
+        ans = 288.0_sp
+        print *, "  dx=0.5", val, ans
+        call check(abs(val - ans) < tol_sp)
+
+        x = [(0.25_sp*(i-1), i = 1, n)]
+        val = simps(y, x)
+        ans = 144.0_sp
+        print *, "  x=0,0.25,0.5,...", val, ans
+        call check(abs(val - ans) < tol_sp)
+    end subroutine test_simps_sp
+
+
+    subroutine test_simps_dp
+        integer, parameter :: n = 13
+        real(dp), dimension(n) :: y
+        real(dp), dimension(n) :: x
+        real(dp) :: val
+        real(dp) :: ans
+        integer :: i
+
+        print *, "test_simps_dp"
+
+        y = [(real(i-1, dp)**2, i = 1, n)]
+
+        val = simps(y, 1.0_dp)
+        ans = 576.0_dp
+        print *, "  dx=1", val, ans
+        call check(abs(val - ans) < tol_dp)
+
+        val = simps(y, 0.5_dp)
+        ans = 288.0_dp
+        print *, "  dx=0.5", val, ans
+        call check(abs(val - ans) < tol_dp)
+
+        x = [(0.25_dp*(i-1), i = 1, n)]
+        val = simps(y, x)
+        ans = 144.0_dp
+        print *, "  x=0,0.25,0.5,...", val, ans
+        call check(abs(val - ans) < tol_dp)
+    end subroutine test_simps_dp
+
+
+    subroutine test_simps_qp
+        integer, parameter :: n = 13
+        real(qp), dimension(n) :: y
+        real(qp), dimension(n) :: x
+        real(qp) :: val
+        real(qp) :: ans
+        integer :: i
+
+        print *, "test_simps_qp"
+
+        y = [(real(i-1, qp)**2, i = 1, n)]
+
+        val = simps(y, 1.0_qp)
+        ans = 576.0_qp
+        print *, "  dx=1", val, ans
+        call check(abs(val - ans) < tol_qp)
+
+        val = simps(y, 0.5_qp)
+        ans = 288.0_qp
+        print *, "  dx=0.5", val, ans
+        call check(abs(val - ans) < tol_qp)
+
+        x = [(0.25_qp*(i-1), i = 1, n)]
+        val = simps(y, x)
+        ans = 144.0_qp
+        print *, "  x=0,0.25,0.5,...", val, ans
+        call check(abs(val - ans) < tol_qp)
+    end subroutine test_simps_qp
+
+
+    subroutine test_simps_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_simps_weights_sp"
+
+        y = [(real(i-1, sp), i = 1, n)]
+
+        x = y
+        w = simps_weights(x)
+        val = sum(w*y)
+        ans = simps(y, x)
+        print *, "  ", val, ans
+        call check(abs(val - ans) < tol_sp)
+    end subroutine test_simps_weights_sp
+
+
+    subroutine test_simps_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_simps_weights_dp"
+
+        y = [(real(i-1, dp), i = 1, n)]
+
+        x = y
+        w = simps_weights(x)
+        val = sum(w*y)
+        ans = simps(y, x)
+        print *, "  ", val, ans
+        call check(abs(val - ans) < tol_dp)
+    end subroutine test_simps_weights_dp
+
+
+    subroutine test_simps_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_simps_weights_qp"
+
+        y = [(real(i-1, qp), i = 1, n)]
+
+        x = y
+        w = simps_weights(x)
+        val = sum(w*y)
+        ans = simps(y, x)
+        print *, "  ", val, ans
+        call check(abs(val - ans) < tol_qp)
+    end subroutine test_simps_weights_qp
+
+
+    subroutine test_simps_zero_sp
+        real(sp), dimension(0) :: a
+
+        print *, "test_simps_zero_sp"
+
+        call check(abs(simps(a, 1.0_sp)) < epsilon(0.0_sp))
+        call check(abs(simps([1.0_sp], 1.0_sp)) < epsilon(0.0_sp))
+        call check(abs(simps(a, a)) < epsilon(0.0_sp))
+        call check(abs(simps([1.0_sp], [1.0_sp])) < epsilon(0.0_sp))
+    end subroutine test_simps_zero_sp
+
+
+    subroutine test_simps_zero_dp
+        real(dp), dimension(0) :: a
+
+        print *, "test_simps_zero_dp"
+
+        call check(abs(simps(a, 1.0_dp)) < epsilon(0.0_dp))
+        call check(abs(simps([1.0_dp], 1.0_dp)) < epsilon(0.0_dp))
+        call check(abs(simps(a, a)) < epsilon(0.0_dp))
+        call check(abs(simps([1.0_dp], [1.0_dp])) < epsilon(0.0_dp))
+    end subroutine test_simps_zero_dp
+
+
+    subroutine test_simps_zero_qp
+        real(qp), dimension(0) :: a
+
+        print *, "test_simps_zero_qp"
+
+        call check(abs(simps(a, 1.0_qp)) < epsilon(0.0_qp))
+        call check(abs(simps([1.0_qp], 1.0_qp)) < epsilon(0.0_qp))
+        call check(abs(simps(a, a)) < epsilon(0.0_qp))
+        call check(abs(simps([1.0_qp], [1.0_qp])) < epsilon(0.0_qp))
+    end subroutine test_simps_zero_qp
+
+
+    subroutine test_simps_even_sp
+        integer, parameter :: n = 11
+        real(sp), dimension(n) :: y
+        real(sp), dimension(n) :: x
+        real(sp) :: val
+        real(sp) :: ans
+        integer :: i
+        integer :: even
+
+        print *, "test_simps_even_sp"
+
+        y = [(3.0_sp*real(i-1, sp)**2, i = 1, n)]
+
+        do even = -1, 1
+            print *, "even=", even
+
+            val = simps(y, 1.0_sp)
+            ans = 1000.0_sp
+            print *, "  dx=1", val, ans
+            call check(abs(val - ans) < tol_sp)
+
+            val = simps(y, 0.5_sp)
+            ans = 500.0_sp
+            print *, "  dx=0.5", val, ans
+            call check(abs(val - ans) < tol_sp)
+
+            x = [(0.25_sp*(i-1), i = 1, n)]
+            val = simps(y, x)
+            ans = 250.0_sp
+            print *, "  x=0,0.25,0.5,...", val, ans
+            call check(abs(val - ans) < tol_sp)
+        end do
+    end subroutine test_simps_even_sp
+
+
+    subroutine test_simps_even_dp
+        integer, parameter :: n = 11
+        real(dp), dimension(n) :: y
+        real(dp), dimension(n) :: x
+        real(dp) :: val
+        real(dp) :: ans
+        integer :: i
+
+        print *, "test_simps_even_dp"
+
+        y = [(3.0_dp*real(i-1, dp)**2, i = 1, n)]
+
+        val = simps(y, 1.0_dp)
+        ans = 1000.0_dp
+        print *, "  dx=1", val, ans
+        call check(abs(val - ans) < tol_dp)
+
+        val = simps(y, 0.5_dp)
+        ans = 500.0_dp
+        print *, "  dx=0.5", val, ans
+        call check(abs(val - ans) < tol_dp)
+
+        x = [(0.25_dp*(i-1), i = 1, n)]
+        val = simps(y, x)
+        ans = 250.0_dp
+        print *, "  x=0,0.25,0.5,...", val, ans
+        call check(abs(val - ans) < tol_dp)
+    end subroutine test_simps_even_dp
+
+
+    subroutine test_simps_even_qp
+        integer, parameter :: n = 11
+        real(qp), dimension(n) :: y
+        real(qp), dimension(n) :: x
+        real(qp) :: val
+        real(qp) :: ans
+        integer :: i
+        integer :: even
+
+        print *, "test_simps_even_qp"
+
+        y = [(3.0_qp*real(i-1, qp)**2, i = 1, n)]
+
+        do even = -1, 1
+            print *, "  even=", even
+
+            val = simps(y, 1.0_qp)
+            ans = 1000.0_qp
+            print *, "  dx=1", val, ans
+            call check(abs(val - ans) < tol_qp)
+
+            val = simps(y, 0.5_qp)
+            ans = 500.0_qp
+            print *, "  dx=0.5", val, ans
+            call check(abs(val - ans) < tol_qp)
+
+            x = [(0.25_qp*(i-1), i = 1, n)]
+            val = simps(y, x)
+            ans = 250.0_qp
+            print *, "  x=0,0.25,0.5,...", val, ans
+            call check(abs(val - ans) < tol_qp)
+        end do
+    end subroutine test_simps_even_qp
+
+
+    subroutine test_simps_weights_even_sp
+        integer, parameter :: n = 16
+        real(sp), dimension(n) :: y
+        real(sp), dimension(n) :: x
+        real(sp), dimension(n) :: w
+        integer :: i
+        real(sp) :: val
+        real(sp) :: ans
+        integer :: even
+
+        print *, "test_simps_weights_even_sp"
+
+        y = [(real(i-1, sp), i = 1, n)]
+        x = y
+
+        do even = -1, 1
+            w = simps_weights(x)
+            val = sum(w*y)
+            ans = simps(y, x)
+            print *, "  even=", even, val, ans
+            call check(abs(val - ans) < tol_sp)
+        end do
+    end subroutine test_simps_weights_even_sp
+
+
+    subroutine test_simps_weights_even_dp
+        integer, parameter :: n = 16
+        real(dp), dimension(n) :: y
+        real(dp), dimension(n) :: x
+        real(dp), dimension(n) :: w
+        integer :: i
+        real(dp) :: val
+        real(dp) :: ans
+        integer :: even
+
+        print *, "test_simps_weights_even_dp"
+
+        y = [(real(i-1, dp), i = 1, n)]
+        x = y
+
+        do even = -1, 1
+            w = simps_weights(x)
+            val = sum(w*y)
+            ans = simps(y, x)
+            print *, "  even=", even, val, ans
+            call check(abs(val - ans) < tol_dp)
+        end do
+    end subroutine test_simps_weights_even_dp
+
+
+    subroutine test_simps_weights_even_qp
+        integer, parameter :: n = 16
+        real(qp), dimension(n) :: y
+        real(qp), dimension(n) :: x
+        real(qp), dimension(n) :: w
+        integer :: i
+        real(qp) :: val
+        real(qp) :: ans
+        integer :: even
+
+        print *, "test_simps_weights_even_qp"
+
+        y = [(real(i-1, qp), i = 1, n)]
+
+        x = y
+
+        do even = -1, 1
+            w = simps_weights(x)
+            val = sum(w*y)
+            ans = simps(y, x)
+            print *, "  even=", even, val, ans
+            call check(abs(val - ans) < tol_qp)
+        end do
+    end subroutine test_simps_weights_even_qp
+
+
+    subroutine test_simps_six_sp
+        integer, parameter :: n = 6
+        real(sp), dimension(n) :: y
+        real(sp), dimension(n) :: x
+        real(sp) :: val
+        real(sp) :: ans
+        integer :: i
+        integer :: even
+
+        print *, "test_simps_six_sp"
+
+        y = [(3.0_sp*real(i-1, sp)**2, i = 1, n)]
+
+        do even = -1, 1
+            print *, "even=", even
+
+            val = simps(y, 1.0_sp)
+            ans = 125.0_sp
+            print *, "  dx=1", val, ans
+            call check(abs(val - ans) < tol_sp)
+
+            val = simps(y, 0.5_sp)
+            ans = 62.5_sp
+            print *, "  dx=0.5", val, ans
+            call check(abs(val - ans) < tol_sp)
+
+            x = [(0.25_sp*(i-1), i = 1, n)]
+            val = simps(y, x)
+            ans = 31.25_sp
+            print *, "  x=0,0.25,0.5,...", val, ans
+            call check(abs(val - ans) < tol_sp)
+        end do
+    end subroutine test_simps_six_sp
+
+
+    subroutine test_simps_six_dp
+        integer, parameter :: n = 6
+        real(dp), dimension(n) :: y
+        real(dp), dimension(n) :: x
+        real(dp) :: val
+        real(dp) :: ans
+        integer :: i
+
+        print *, "test_simps_six_dp"
+
+        y = [(3.0_dp*real(i-1, dp)**2, i = 1, n)]
+
+        val = simps(y, 1.0_dp)
+        ans = 125.0_dp
+        print *, "  dx=1", val, ans
+        call check(abs(val - ans) < tol_dp)
+
+        val = simps(y, 0.5_dp)
+        ans = 62.5_dp
+        print *, "  dx=0.5", val, ans
+        call check(abs(val - ans) < tol_dp)
+
+        x = [(0.25_dp*(i-1), i = 1, n)]
+        val = simps(y, x)
+        ans = 31.25_dp
+        print *, "  x=0,0.25,0.5,...", val, ans
+        call check(abs(val - ans) < tol_dp)
+    end subroutine test_simps_six_dp
+
+
+    subroutine test_simps_six_qp
+        integer, parameter :: n = 6
+        real(qp), dimension(n) :: y
+        real(qp), dimension(n) :: x
+        real(qp) :: val
+        real(qp) :: ans
+        integer :: i
+        integer :: even
+
+        print *, "test_simps_six_qp"
+
+        y = [(3.0_qp*real(i-1, qp)**2, i = 1, n)]
+
+        do even = -1, 1
+            print *, "  even=", even
+
+            val = simps(y, 1.0_qp)
+            ans = 125.0_qp
+            print *, "  dx=1", val, ans
+            call check(abs(val - ans) < tol_qp)
+
+            val = simps(y, 0.5_qp)
+            ans = 62.5_qp
+            print *, "  dx=0.5", val, ans
+            call check(abs(val - ans) < tol_qp)
+
+            x = [(0.25_qp*(i-1), i = 1, n)]
+            val = simps(y, x)
+            ans = 31.25_qp
+            print *, "  x=0,0.25,0.5,...", val, ans
+            call check(abs(val - ans) < tol_qp)
+        end do
+    end subroutine test_simps_six_qp
+
+end program