From 56822b586cd7514c4d8386fa13e72c88a77db03e Mon Sep 17 00:00:00 2001
From: "Harris M. Snyder" <harris.snyder@gmail.com>
Date: Sun, 7 Feb 2021 23:21:08 -0500
Subject: [PATCH 01/18] - Added preliminary implementation for Legendre
 polynomials - Added preliminary implementation for Legendre polynomial first
 derivatives - Added preliminary implementation for Gauss-Legendre quadrature
 points/weights - Added preliminary implementation for Gauss-Legendre-Lobatto
 quadrature points/weights

---
 src/CMakeLists.txt                |   3 +
 src/Makefile.manual               |   6 ++
 src/stdlib_functions.f90          |  34 +++++++++
 src/stdlib_functions_legendre.f90 |  73 +++++++++++++++++++
 src/stdlib_quadrature.fypp        |  26 +++++++
 src/stdlib_quadrature_gauss.f90   | 115 ++++++++++++++++++++++++++++++
 6 files changed, 257 insertions(+)
 create mode 100644 src/stdlib_functions.f90
 create mode 100644 src/stdlib_functions_legendre.f90
 create mode 100644 src/stdlib_quadrature_gauss.f90

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 02413b415..4546c5058 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -41,6 +41,9 @@ set(SRC
     stdlib_kinds.f90
     stdlib_logger.f90
     stdlib_system.F90
+    stdlib_functions.f90
+    stdlib_functions_legendre.f90
+    stdlib_quadrature_gauss.f90
     ${outFiles}
 )
 
diff --git a/src/Makefile.manual b/src/Makefile.manual
index dd6f12708..92ffe9a76 100644
--- a/src/Makefile.manual
+++ b/src/Makefile.manual
@@ -4,6 +4,8 @@ SRC = f18estop.f90 \
       stdlib_bitsets_64.f90 \
       stdlib_bitsets_large.f90 \
       stdlib_error.f90 \
+      stdlib_functions.f90 \
+      stdlib_functions_legendre.f90 \
       stdlib_io.f90 \
       stdlib_kinds.f90 \
       stdlib_linalg.f90 \
@@ -12,6 +14,7 @@ SRC = f18estop.f90 \
       stdlib_optval.f90 \
       stdlib_quadrature.f90 \
       stdlib_quadrature_trapz.f90 \
+      stdlib_quadrature_gauss.f90 \
       stdlib_stats.f90 \
       stdlib_stats_mean.f90 \
       stdlib_stats_moment.f90 \
@@ -50,6 +53,8 @@ stdlib_bitsets.o: stdlib_kinds.o
 stdlib_bitsets_64.o: stdlib_bitsets.o
 stdlib_bitsets_large.o: stdlib_bitsets.o
 stdlib_error.o: stdlib_optval.o
+stdlib_functions.o: stdlib_kinds.o
+stdlib_functions_legendre.o: stdlib_kinds.o
 stdlib_io.o: \
 	stdlib_error.o \
 	stdlib_optval.o \
@@ -58,6 +63,7 @@ stdlib_linalg_diag.o: stdlib_kinds.o
 stdlib_logger.o: stdlib_ascii.o stdlib_optval.o
 stdlib_optval.o: stdlib_kinds.o
 stdlib_quadrature.o: stdlib_kinds.o
+stdlib_quadrature_gauss.o: stdlib_kinds.o
 stdlib_stats_mean.o: \
 	stdlib_optval.o \
 	stdlib_kinds.o \
diff --git a/src/stdlib_functions.f90 b/src/stdlib_functions.f90
new file mode 100644
index 000000000..46b6a0dc6
--- /dev/null
+++ b/src/stdlib_functions.f90
@@ -0,0 +1,34 @@
+module stdlib_functions
+    use stdlib_kinds, only: sp, dp, qp
+
+    implicit none
+
+    private
+
+    public :: legendre 
+    public :: dlegendre 
+
+
+    interface legendre
+        !! version: experimental
+        !! 
+        !! Legendre polynomial
+        pure elemental module function legendre_fp64(N,x) result(L)
+            integer, intent(in) :: N
+            real(dp), intent(in) :: x
+            real(dp) :: L
+        end function
+    end interface
+
+    interface dlegendre
+        !! version: experimental
+        !! 
+        !! First derivative Legendre polynomial
+        pure elemental module function dlegendre_fp64(N,x) result(DL)
+            integer, intent(in) :: N
+            real(dp), intent(in) :: x
+            real(dp) :: DL
+        end function
+    end interface
+
+end module stdlib_functions
diff --git a/src/stdlib_functions_legendre.f90 b/src/stdlib_functions_legendre.f90
new file mode 100644
index 000000000..b9532c95a
--- /dev/null
+++ b/src/stdlib_functions_legendre.f90
@@ -0,0 +1,73 @@
+submodule (stdlib_functions) stdlib_functions_legendre
+    use stdlib_kinds, only: sp, dp, qp
+    implicit none
+
+contains
+
+    ! derivatives of Legendre polynomials
+    ! unspecified behaviour if N is negative
+    pure elemental module function dlegendre_fp64(N,x) result(DL)
+        integer, intent(in) :: N
+        real(dp), intent(in) :: x
+        real(dp) :: DL
+
+        select case(N)
+            case(0)
+                DL = 0
+            case(1)
+                DL = 1
+            case default
+                block
+                    real(dp) :: L_down1, L_down2, L
+                    real(dp) :: DL_down1, DL_down2
+                    integer :: i 
+
+                    L_down1  = x
+                    DL_down1 = 1
+
+                    L_down2  = 1
+                    DL_down2 = 0
+
+                    do i = 2, N
+                        L = (2*i-1)*x*L_down1/i - (i-1)*L_down2/i
+                        DL = DL_down2 + (2*i-1)*L_down1
+                        L_down2 = L_down1
+                        L_down1 = L
+                        DL_down2 = DL_down1
+                        DL_down1 = DL
+                    end do
+                end block
+        end select
+    end function
+
+    ! Legendre polynomials
+    ! unspecified behaviour if N is negative
+    pure elemental module function legendre_fp64(N,x) result(L)
+        integer, intent(in) :: N
+        real(dp), intent(in) :: x
+        real(dp) :: L
+        select case(N)
+            case(0)
+                L  = 1
+            case(1)
+                L  = x
+            case default
+                block
+                    real(dp) :: L_down1, L_down2
+                    integer :: i 
+
+                    L_down1  = x
+                    L_down2  = 1
+
+                    do i = 2, N
+                        L = (2*i-1)*x*L_down1/i - (i-1)*L_down2/i
+                        L_down2 = L_down1
+                        L_down1 = L
+                    end do
+                end block
+        end select
+    end function
+
+end submodule
+
+
diff --git a/src/stdlib_quadrature.fypp b/src/stdlib_quadrature.fypp
index e8c7e0531..63d99c027 100644
--- a/src/stdlib_quadrature.fypp
+++ b/src/stdlib_quadrature.fypp
@@ -12,6 +12,8 @@ module stdlib_quadrature
     public :: trapz_weights
     public :: simps
     public :: simps_weights
+    public :: gauss_legendre
+    public :: gauss_legendre_lobatto
 
 
     interface trapz
@@ -90,6 +92,28 @@ module stdlib_quadrature
     end interface simps_weights
 
 
+    interface gauss_legendre
+        !! version: experimental
+        !!
+        !! Computes Gauss-Legendre quadrature nodes and weights.
+        pure module subroutine gauss_legendre_fp64 (x, w, interval)
+            real(dp), intent(out) :: x(:), w(:)
+            real(dp), intent(in), optional :: interval(2)
+        end subroutine
+    end interface gauss_legendre
+
+
+    interface gauss_legendre_lobatto
+        !! version: experimental
+        !!
+        !! Computes Gauss-Legendre-Lobatto quadrature nodes and weights.
+        pure module subroutine gauss_legendre_lobatto_fp64 (x, w, interval)
+            real(dp), intent(out) :: x(:), w(:)
+            real(dp), intent(in), optional :: interval(2)
+        end subroutine
+    end interface gauss_legendre_lobatto
+
+
     ! 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.
@@ -103,4 +127,6 @@ module stdlib_quadrature
         #:endfor
     end interface
 
+
+
 end module stdlib_quadrature
diff --git a/src/stdlib_quadrature_gauss.f90 b/src/stdlib_quadrature_gauss.f90
new file mode 100644
index 000000000..e80b8fbee
--- /dev/null
+++ b/src/stdlib_quadrature_gauss.f90
@@ -0,0 +1,115 @@
+submodule (stdlib_quadrature) stdlib_quadrature_gauss
+    use stdlib_kinds, only: sp, dp, qp
+    use stdlib_functions, only: legendre, dlegendre
+    implicit none
+
+    real(dp), parameter :: PI = acos(-1._dp)
+    real(dp), parameter :: tolerance = 4._dp * epsilon(1._dp)
+    integer, parameter :: newton_iters = 100
+
+contains
+
+    pure module subroutine gauss_legendre_fp64 (x, w, interval)
+        real(dp), intent(out) :: x(:), w(:)
+        real(dp), intent(in), optional :: interval(2)
+
+        associate (N => size(x)-1 )
+        select case (N)
+            case (0)
+                x = 0._dp
+                w = 2._dp
+            case (1)
+                x = [-sqrt(1._dp/3._dp), sqrt(1._dp/3._dp)]
+                w = [1._dp, 1._dp]
+            case default
+                block
+                integer :: i,j
+                real(dp) :: leg, dleg, delta
+
+                do i = 0, int(floor((N+1)/2._dp)-1)
+                    x(i+1) = -cos((2*i+1)/(2._dp*N+2._dp) * PI)
+                    do j = 0, newton_iters-1
+                        leg  = legendre(N+1,x(i+1))
+                        dleg = dlegendre(N+1,x(i+1))
+                        delta = -leg/dleg
+                        x(i+1) = x(i+1) + delta
+                        if ( abs(delta) <= tolerance * abs(x(i+1)) )  exit
+                    end do
+                    x(N-i+1) = -x(i+1)
+
+                    dleg = dlegendre(N+1,x(i+1))
+                    w(i+1)   = 2._dp/((1-x(i+1)**2)*dleg**2) 
+                    w(N-i+1) = w(i+1)
+                end do
+
+                if (mod(N,2) == 0) then
+                    x(N/2+1) = 0.0
+
+                    dleg = dlegendre(N+1, 0.0_dp)
+                    w(N/2+1) = 2._dp/(dleg**2) 
+                end if
+                end block
+        end select
+        end associate
+
+        if (present(interval)) then
+            associate ( a => interval(1) , b => interval(2) )
+            x = 0.5*(b-a)*x+0.5*(b+a)
+            w = 0.5*(b-a)*w
+            end associate
+        end if
+    end subroutine
+
+    pure module subroutine gauss_legendre_lobatto_fp64 (x, w, interval)
+        real(dp), intent(out) :: x(:), w(:)
+        real(dp), intent(in), optional :: interval(2)
+
+        associate (N => size(x)-1)
+        select case (N)
+            case (1)
+                x = [-1._dp, 1._dp]
+                w = [ 1._dp, 1._dp]
+            case default
+                block
+                integer :: i,j
+                real(dp) :: leg, dleg, delta
+
+                x(1)   = -1._dp
+                x(N+1) =  1._dp
+                w(1)   =  2._dp/(N*(N+1._dp))
+                w(N+1) =  2._dp/(N*(N+1._dp))
+
+                do i = 1, int(floor((N+1)/2._dp)-1)
+                    x(i+1) = -cos( (i+0.25_dp)*PI/N  - 3/(8*N*PI*(i+0.25_dp)))
+                    do j = 0, newton_iters-1
+                        leg  = legendre(N+1,x(i+1)) - legendre(N-1,x(i+1))
+                        dleg = dlegendre(N+1,x(i+1)) - dlegendre(N-1,x(i+1))
+                        delta = -leg/dleg
+                        x(i+1) = x(i+1) + delta
+                        if ( abs(delta) <= tolerance * abs(x(i+1)) )  exit
+                    end do
+                    x(N-i+1) = -x(i+1)
+
+                    leg = legendre(N, x(i+1))
+                    w(i+1)   = 2._dp/(N*(N+1._dp)*leg**2) 
+                    w(N-i+1) = w(i+1)
+                end do
+
+                if (mod(N,2) == 0) then
+                    x(N/2+1) = 0.0
+
+                    leg = legendre(N, 0.0_dp)
+                    w(N/2+1)   = 2._dp/(N*(N+1._dp)*leg**2) 
+                end if
+                end block
+        end select
+        end associate
+        
+        if (present(interval)) then
+            associate ( a => interval(1) , b => interval(2) )
+            x = 0.5*(b-a)*x+0.5*(b+a)
+            w = 0.5*(b-a)*w
+            end associate
+        end if
+    end subroutine
+end submodule    

From 1a914ddbbb5956c696a8098d48a19712f27d5ad9 Mon Sep 17 00:00:00 2001
From: Harris Snyder <harris.snyder@gmail.com>
Date: Mon, 8 Feb 2021 13:28:39 -0500
Subject: [PATCH 02/18] Update src/stdlib_quadrature_gauss.f90

Co-authored-by: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com>
---
 src/stdlib_quadrature_gauss.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/stdlib_quadrature_gauss.f90 b/src/stdlib_quadrature_gauss.f90
index e80b8fbee..711067462 100644
--- a/src/stdlib_quadrature_gauss.f90
+++ b/src/stdlib_quadrature_gauss.f90
@@ -26,7 +26,7 @@ pure module subroutine gauss_legendre_fp64 (x, w, interval)
                 integer :: i,j
                 real(dp) :: leg, dleg, delta
 
-                do i = 0, int(floor((N+1)/2._dp)-1)
+                do i = 0, (N+1)/2-1
                     x(i+1) = -cos((2*i+1)/(2._dp*N+2._dp) * PI)
                     do j = 0, newton_iters-1
                         leg  = legendre(N+1,x(i+1))

From 5937f90595e41fbcfd7c2300b83cc1fc0cd1d21d Mon Sep 17 00:00:00 2001
From: "Harris M. Snyder" <harris.snyder@gmail.com>
Date: Mon, 8 Feb 2021 14:18:56 -0500
Subject: [PATCH 03/18] Some changes in accorance with preliminary review
 comments

---
 src/stdlib_functions.f90          | 12 ++---
 src/stdlib_functions_legendre.f90 | 72 ++++++++++++-------------
 src/stdlib_quadrature_gauss.f90   | 88 +++++++++++++++++--------------
 3 files changed, 89 insertions(+), 83 deletions(-)

diff --git a/src/stdlib_functions.f90 b/src/stdlib_functions.f90
index 46b6a0dc6..6e36cd789 100644
--- a/src/stdlib_functions.f90
+++ b/src/stdlib_functions.f90
@@ -13,10 +13,10 @@ module stdlib_functions
         !! version: experimental
         !! 
         !! Legendre polynomial
-        pure elemental module function legendre_fp64(N,x) result(L)
-            integer, intent(in) :: N
+        pure elemental module function legendre_fp64(n,x) result(leg)
+            integer, intent(in) :: n
             real(dp), intent(in) :: x
-            real(dp) :: L
+            real(dp) :: leg
         end function
     end interface
 
@@ -24,10 +24,10 @@ pure elemental module function legendre_fp64(N,x) result(L)
         !! version: experimental
         !! 
         !! First derivative Legendre polynomial
-        pure elemental module function dlegendre_fp64(N,x) result(DL)
-            integer, intent(in) :: N
+        pure elemental module function dlegendre_fp64(n,x) result(dleg)
+            integer, intent(in) :: n
             real(dp), intent(in) :: x
-            real(dp) :: DL
+            real(dp) :: dleg
         end function
     end interface
 
diff --git a/src/stdlib_functions_legendre.f90 b/src/stdlib_functions_legendre.f90
index b9532c95a..195de0b16 100644
--- a/src/stdlib_functions_legendre.f90
+++ b/src/stdlib_functions_legendre.f90
@@ -4,65 +4,65 @@
 
 contains
 
-    ! derivatives of Legendre polynomials
-    ! unspecified behaviour if N is negative
-    pure elemental module function dlegendre_fp64(N,x) result(DL)
-        integer, intent(in) :: N
+    ! derivatives of legegendre polynomials
+    ! unspecified behaviour if n is negative
+    pure elemental module function dlegendre_fp64(n,x) result(dleg)
+        integer, intent(in) :: n
         real(dp), intent(in) :: x
-        real(dp) :: DL
+        real(dp) :: dleg
 
-        select case(N)
+        select case(n)
             case(0)
-                DL = 0
+                dleg = 0
             case(1)
-                DL = 1
+                dleg = 1
             case default
                 block
-                    real(dp) :: L_down1, L_down2, L
-                    real(dp) :: DL_down1, DL_down2
+                    real(dp) :: leg_down1, leg_down2, leg
+                    real(dp) :: dleg_down1, dleg_down2
                     integer :: i 
 
-                    L_down1  = x
-                    DL_down1 = 1
+                    leg_down1  = x
+                    dleg_down1 = 1
 
-                    L_down2  = 1
-                    DL_down2 = 0
+                    leg_down2  = 1
+                    dleg_down2 = 0
 
-                    do i = 2, N
-                        L = (2*i-1)*x*L_down1/i - (i-1)*L_down2/i
-                        DL = DL_down2 + (2*i-1)*L_down1
-                        L_down2 = L_down1
-                        L_down1 = L
-                        DL_down2 = DL_down1
-                        DL_down1 = DL
+                    do i = 2, n
+                        leg = (2*i-1)*x*leg_down1/i - (i-1)*leg_down2/i
+                        dleg = dleg_down2 + (2*i-1)*leg_down1
+                        leg_down2 = leg_down1
+                        leg_down1 = leg
+                        dleg_down2 = dleg_down1
+                        dleg_down1 = dleg
                     end do
                 end block
         end select
     end function
 
-    ! Legendre polynomials
-    ! unspecified behaviour if N is negative
-    pure elemental module function legendre_fp64(N,x) result(L)
-        integer, intent(in) :: N
+    ! legegendre polynomials
+    ! unspecified behaviour if n is negative
+    pure elemental module function legendre_fp64(n,x) result(leg)
+        integer, intent(in) :: n
         real(dp), intent(in) :: x
-        real(dp) :: L
-        select case(N)
+        real(dp) :: leg
+        select case(n)
             case(0)
-                L  = 1
+                leg  = 1
             case(1)
-                L  = x
+                leg  = x
             case default
                 block
-                    real(dp) :: L_down1, L_down2
+                    real(dp) :: leg_down1, leg_down2
                     integer :: i 
 
-                    L_down1  = x
-                    L_down2  = 1
+                    leg_down1  = x
+                    leg_down2  = 1
 
-                    do i = 2, N
-                        L = (2*i-1)*x*L_down1/i - (i-1)*L_down2/i
-                        L_down2 = L_down1
-                        L_down1 = L
+                    do i = 2, n
+                        leg = (2*i-1)*x*leg_down1/i - (i-1)*leg_down2/i
+                        leg_down2 = leg_down1
+                        leg_down1 = leg
                     end do
                 end block
         end select
diff --git a/src/stdlib_quadrature_gauss.f90 b/src/stdlib_quadrature_gauss.f90
index e80b8fbee..7a3c80c65 100644
--- a/src/stdlib_quadrature_gauss.f90
+++ b/src/stdlib_quadrature_gauss.f90
@@ -3,7 +3,7 @@
     use stdlib_functions, only: legendre, dlegendre
     implicit none
 
-    real(dp), parameter :: PI = acos(-1._dp)
+    real(dp), parameter :: pi = acos(-1._dp)
     real(dp), parameter :: tolerance = 4._dp * epsilon(1._dp)
     integer, parameter :: newton_iters = 100
 
@@ -13,40 +13,41 @@ pure module subroutine gauss_legendre_fp64 (x, w, interval)
         real(dp), intent(out) :: x(:), w(:)
         real(dp), intent(in), optional :: interval(2)
 
-        associate (N => size(x)-1 )
-        select case (N)
+        associate (n => size(x)-1 )
+        select case (n)
             case (0)
-                x = 0._dp
-                w = 2._dp
+                x = 0
+                w = 2
             case (1)
-                x = [-sqrt(1._dp/3._dp), sqrt(1._dp/3._dp)]
-                w = [1._dp, 1._dp]
+                x(1) = -sqrt(1._dp/3._dp)
+                x(2) = -x(1)
+                w = 1
             case default
                 block
                 integer :: i,j
                 real(dp) :: leg, dleg, delta
 
-                do i = 0, int(floor((N+1)/2._dp)-1)
-                    x(i+1) = -cos((2*i+1)/(2._dp*N+2._dp) * PI)
+                do i = 0, (n+1)/2 - 1
+                    x(i+1) = -cos((2*i+1)/(2._dp*n+2._dp) * pi)
                     do j = 0, newton_iters-1
-                        leg  = legendre(N+1,x(i+1))
-                        dleg = dlegendre(N+1,x(i+1))
+                        leg  = legendre(n+1,x(i+1))
+                        dleg = dlegendre(n+1,x(i+1))
                         delta = -leg/dleg
                         x(i+1) = x(i+1) + delta
                         if ( abs(delta) <= tolerance * abs(x(i+1)) )  exit
                     end do
-                    x(N-i+1) = -x(i+1)
+                    x(n-i+1) = -x(i+1)
 
-                    dleg = dlegendre(N+1,x(i+1))
+                    dleg = dlegendre(n+1,x(i+1))
                     w(i+1)   = 2._dp/((1-x(i+1)**2)*dleg**2) 
-                    w(N-i+1) = w(i+1)
+                    w(n-i+1) = w(i+1)
                 end do
 
-                if (mod(N,2) == 0) then
-                    x(N/2+1) = 0.0
+                if (mod(n,2) == 0) then
+                    x(n/2+1) = 0
 
-                    dleg = dlegendre(N+1, 0.0_dp)
-                    w(N/2+1) = 2._dp/(dleg**2) 
+                    dleg = dlegendre(n+1, 0.0_dp)
+                    w(n/2+1) = 2._dp/(dleg**2) 
                 end if
                 end block
         end select
@@ -54,8 +55,10 @@ pure module subroutine gauss_legendre_fp64 (x, w, interval)
 
         if (present(interval)) then
             associate ( a => interval(1) , b => interval(2) )
-            x = 0.5*(b-a)*x+0.5*(b+a)
-            w = 0.5*(b-a)*w
+            x = 0.5_dp*(b-a)*x+0.5_dp*(b+a)
+            x(1)       = interval(1)
+            x(size(x)) = interval(2)
+            w = 0.5_dp*(b-a)*w
             end associate
         end if
     end subroutine
@@ -64,42 +67,43 @@ pure module subroutine gauss_legendre_lobatto_fp64 (x, w, interval)
         real(dp), intent(out) :: x(:), w(:)
         real(dp), intent(in), optional :: interval(2)
 
-        associate (N => size(x)-1)
-        select case (N)
+        associate (n => size(x)-1)
+        select case (n)
             case (1)
-                x = [-1._dp, 1._dp]
-                w = [ 1._dp, 1._dp]
+                x(1) = -1
+                x(2) =  1
+                w = 1
             case default
                 block
                 integer :: i,j
                 real(dp) :: leg, dleg, delta
 
                 x(1)   = -1._dp
-                x(N+1) =  1._dp
-                w(1)   =  2._dp/(N*(N+1._dp))
-                w(N+1) =  2._dp/(N*(N+1._dp))
+                x(n+1) =  1._dp
+                w(1)   =  2._dp/(n*(n+1._dp))
+                w(n+1) =  2._dp/(n*(n+1._dp))
 
-                do i = 1, int(floor((N+1)/2._dp)-1)
-                    x(i+1) = -cos( (i+0.25_dp)*PI/N  - 3/(8*N*PI*(i+0.25_dp)))
+                do i = 1, (n+1)/2 - 1
+                    x(i+1) = -cos( (i+0.25_dp)*pi/n  - 3/(8*n*pi*(i+0.25_dp)))
                     do j = 0, newton_iters-1
-                        leg  = legendre(N+1,x(i+1)) - legendre(N-1,x(i+1))
-                        dleg = dlegendre(N+1,x(i+1)) - dlegendre(N-1,x(i+1))
+                        leg  = legendre(n+1,x(i+1)) - legendre(n-1,x(i+1))
+                        dleg = dlegendre(n+1,x(i+1)) - dlegendre(n-1,x(i+1))
                         delta = -leg/dleg
                         x(i+1) = x(i+1) + delta
                         if ( abs(delta) <= tolerance * abs(x(i+1)) )  exit
                     end do
-                    x(N-i+1) = -x(i+1)
+                    x(n-i+1) = -x(i+1)
 
-                    leg = legendre(N, x(i+1))
-                    w(i+1)   = 2._dp/(N*(N+1._dp)*leg**2) 
-                    w(N-i+1) = w(i+1)
+                    leg = legendre(n, x(i+1))
+                    w(i+1)   = 2._dp/(n*(n+1._dp)*leg**2) 
+                    w(n-i+1) = w(i+1)
                 end do
 
-                if (mod(N,2) == 0) then
-                    x(N/2+1) = 0.0
+                if (mod(n,2) == 0) then
+                    x(n/2+1) = 0
 
-                    leg = legendre(N, 0.0_dp)
-                    w(N/2+1)   = 2._dp/(N*(N+1._dp)*leg**2) 
+                    leg = legendre(n, 0.0_dp)
+                    w(n/2+1)   = 2._dp/(n*(n+1._dp)*leg**2) 
                 end if
                 end block
         end select
@@ -107,8 +111,10 @@ pure module subroutine gauss_legendre_lobatto_fp64 (x, w, interval)
         
         if (present(interval)) then
             associate ( a => interval(1) , b => interval(2) )
-            x = 0.5*(b-a)*x+0.5*(b+a)
-            w = 0.5*(b-a)*w
+            x = 0.5_dp*(b-a)*x+0.5_dp*(b+a)
+            x(1)       = interval(1)
+            x(size(x)) = interval(2)
+            w = 0.5_dp*(b-a)*w
             end associate
         end if
     end subroutine

From e0add1a260e1ba115733e77a96e8b19c0b847fd9 Mon Sep 17 00:00:00 2001
From: "Harris M. Snyder" <harris.snyder@gmail.com>
Date: Tue, 9 Feb 2021 10:44:09 -0500
Subject: [PATCH 04/18] changed stdlib_functions to stdlib_specialfunctions

---
 src/CMakeLists.txt                                        | 4 ++--
 src/Makefile.manual                                       | 8 ++++----
 src/stdlib_quadrature_gauss.f90                           | 3 ++-
 src/{stdlib_functions.f90 => stdlib_specialfunctions.f90} | 4 ++--
 ..._legendre.f90 => stdlib_specialfunctions_legendre.f90} | 2 +-
 5 files changed, 11 insertions(+), 10 deletions(-)
 rename src/{stdlib_functions.f90 => stdlib_specialfunctions.f90} (92%)
 rename src/{stdlib_functions_legendre.f90 => stdlib_specialfunctions_legendre.f90} (96%)

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 4546c5058..33d103e40 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -41,8 +41,8 @@ set(SRC
     stdlib_kinds.f90
     stdlib_logger.f90
     stdlib_system.F90
-    stdlib_functions.f90
-    stdlib_functions_legendre.f90
+    stdlib_specialfunctions.f90
+    stdlib_specialfunctions_legendre.f90
     stdlib_quadrature_gauss.f90
     ${outFiles}
 )
diff --git a/src/Makefile.manual b/src/Makefile.manual
index 92ffe9a76..8a0572713 100644
--- a/src/Makefile.manual
+++ b/src/Makefile.manual
@@ -4,8 +4,8 @@ SRC = f18estop.f90 \
       stdlib_bitsets_64.f90 \
       stdlib_bitsets_large.f90 \
       stdlib_error.f90 \
-      stdlib_functions.f90 \
-      stdlib_functions_legendre.f90 \
+      stdlib_specialfunctions.f90 \
+      stdlib_specialfunctions_legendre.f90 \
       stdlib_io.f90 \
       stdlib_kinds.f90 \
       stdlib_linalg.f90 \
@@ -53,8 +53,8 @@ stdlib_bitsets.o: stdlib_kinds.o
 stdlib_bitsets_64.o: stdlib_bitsets.o
 stdlib_bitsets_large.o: stdlib_bitsets.o
 stdlib_error.o: stdlib_optval.o
-stdlib_functions.o: stdlib_kinds.o
-stdlib_functions_legendre.o: stdlib_kinds.o
+stdlib_specialfunctions.o: stdlib_kinds.o
+stdlib_specialfunctions_legendre.o: stdlib_kinds.o
 stdlib_io.o: \
 	stdlib_error.o \
 	stdlib_optval.o \
diff --git a/src/stdlib_quadrature_gauss.f90 b/src/stdlib_quadrature_gauss.f90
index 7a3c80c65..8f909f3ea 100644
--- a/src/stdlib_quadrature_gauss.f90
+++ b/src/stdlib_quadrature_gauss.f90
@@ -1,6 +1,6 @@
 submodule (stdlib_quadrature) stdlib_quadrature_gauss
     use stdlib_kinds, only: sp, dp, qp
-    use stdlib_functions, only: legendre, dlegendre
+    use stdlib_specialfunctions, only: legendre, dlegendre
     implicit none
 
     real(dp), parameter :: pi = acos(-1._dp)
@@ -84,6 +84,7 @@ pure module subroutine gauss_legendre_lobatto_fp64 (x, w, interval)
                 w(n+1) =  2._dp/(n*(n+1._dp))
 
                 do i = 1, (n+1)/2 - 1
+                    ! initial guess from an approximate form given by SV Parter (1999)
                     x(i+1) = -cos( (i+0.25_dp)*pi/n  - 3/(8*n*pi*(i+0.25_dp)))
                     do j = 0, newton_iters-1
                         leg  = legendre(n+1,x(i+1)) - legendre(n-1,x(i+1))
diff --git a/src/stdlib_functions.f90 b/src/stdlib_specialfunctions.f90
similarity index 92%
rename from src/stdlib_functions.f90
rename to src/stdlib_specialfunctions.f90
index 6e36cd789..32859fe6d 100644
--- a/src/stdlib_functions.f90
+++ b/src/stdlib_specialfunctions.f90
@@ -1,4 +1,4 @@
-module stdlib_functions
+module stdlib_specialfunctions
     use stdlib_kinds, only: sp, dp, qp
 
     implicit none
@@ -31,4 +31,4 @@ pure elemental module function dlegendre_fp64(n,x) result(dleg)
         end function
     end interface
 
-end module stdlib_functions
+end module stdlib_specialfunctions
diff --git a/src/stdlib_functions_legendre.f90 b/src/stdlib_specialfunctions_legendre.f90
similarity index 96%
rename from src/stdlib_functions_legendre.f90
rename to src/stdlib_specialfunctions_legendre.f90
index 195de0b16..e18ca283b 100644
--- a/src/stdlib_functions_legendre.f90
+++ b/src/stdlib_specialfunctions_legendre.f90
@@ -1,4 +1,4 @@
-submodule (stdlib_functions) stdlib_functions_legendre
+submodule (stdlib_specialfunctions) stdlib_specialfunctions_legendre
     use stdlib_kinds, only: sp, dp, qp
     implicit none
 

From 29ebe9c7c4b53cd3837fe8e671f5827aef0aee27 Mon Sep 17 00:00:00 2001
From: "Harris M. Snyder" <harris.snyder@gmail.com>
Date: Wed, 10 Feb 2021 10:51:20 -0500
Subject: [PATCH 05/18] added a comment regarding the initial guess for
 gauss-legendre quadrature points

---
 src/stdlib_quadrature_gauss.f90 | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/stdlib_quadrature_gauss.f90 b/src/stdlib_quadrature_gauss.f90
index 8f909f3ea..b67b1de7e 100644
--- a/src/stdlib_quadrature_gauss.f90
+++ b/src/stdlib_quadrature_gauss.f90
@@ -28,6 +28,7 @@ pure module subroutine gauss_legendre_fp64 (x, w, interval)
                 real(dp) :: leg, dleg, delta
 
                 do i = 0, (n+1)/2 - 1
+                    ! use Gauss-Chebyshev points as an initial guess
                     x(i+1) = -cos((2*i+1)/(2._dp*n+2._dp) * pi)
                     do j = 0, newton_iters-1
                         leg  = legendre(n+1,x(i+1))

From 7f91f1dae9d17571c41923ecedc5190fde760985 Mon Sep 17 00:00:00 2001
From: "Harris M. Snyder" <harris.snyder@gmail.com>
Date: Wed, 10 Feb 2021 13:20:52 -0500
Subject: [PATCH 06/18] changed newton iteration loop indexing
 (stdlib_quadrature_gauss.f90) in accordance to a review comment

---
 src/stdlib_quadrature_gauss.f90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/stdlib_quadrature_gauss.f90 b/src/stdlib_quadrature_gauss.f90
index b67b1de7e..682de79c8 100644
--- a/src/stdlib_quadrature_gauss.f90
+++ b/src/stdlib_quadrature_gauss.f90
@@ -30,7 +30,7 @@ pure module subroutine gauss_legendre_fp64 (x, w, interval)
                 do i = 0, (n+1)/2 - 1
                     ! use Gauss-Chebyshev points as an initial guess
                     x(i+1) = -cos((2*i+1)/(2._dp*n+2._dp) * pi)
-                    do j = 0, newton_iters-1
+                    do j = 1, newton_iters
                         leg  = legendre(n+1,x(i+1))
                         dleg = dlegendre(n+1,x(i+1))
                         delta = -leg/dleg
@@ -87,7 +87,7 @@ pure module subroutine gauss_legendre_lobatto_fp64 (x, w, interval)
                 do i = 1, (n+1)/2 - 1
                     ! initial guess from an approximate form given by SV Parter (1999)
                     x(i+1) = -cos( (i+0.25_dp)*pi/n  - 3/(8*n*pi*(i+0.25_dp)))
-                    do j = 0, newton_iters-1
+                    do j = 1, newton_iters
                         leg  = legendre(n+1,x(i+1)) - legendre(n-1,x(i+1))
                         dleg = dlegendre(n+1,x(i+1)) - dlegendre(n-1,x(i+1))
                         delta = -leg/dleg

From 56b942319d344baae4f52a260c2b944e915552f7 Mon Sep 17 00:00:00 2001
From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com>
Date: Sat, 10 Apr 2021 15:17:42 +0200
Subject: [PATCH 07/18] Fix build for GCC 10.2

---
 src/stdlib_quadrature_gauss.f90 | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/stdlib_quadrature_gauss.f90 b/src/stdlib_quadrature_gauss.f90
index 682de79c8..0a346db48 100644
--- a/src/stdlib_quadrature_gauss.f90
+++ b/src/stdlib_quadrature_gauss.f90
@@ -1,5 +1,4 @@
 submodule (stdlib_quadrature) stdlib_quadrature_gauss
-    use stdlib_kinds, only: sp, dp, qp
     use stdlib_specialfunctions, only: legendre, dlegendre
     implicit none
 

From 551ec6dc9375dc8548109c0f6064c34502585b6b Mon Sep 17 00:00:00 2001
From: "Harris M. Snyder" <harris.snyder@gmail.com>
Date: Wed, 23 Jun 2021 16:46:09 -0400
Subject: [PATCH 08/18] added first tests for gaussian quadrature functions

---
 src/tests/quadrature/CMakeLists.txt  |  1 +
 src/tests/quadrature/Makefile.manual |  2 ++
 src/tests/quadrature/test_gauss.f90  | 34 ++++++++++++++++++++++++++++
 3 files changed, 37 insertions(+)
 create mode 100644 src/tests/quadrature/test_gauss.f90

diff --git a/src/tests/quadrature/CMakeLists.txt b/src/tests/quadrature/CMakeLists.txt
index d8069fb86..c8fc618db 100644
--- a/src/tests/quadrature/CMakeLists.txt
+++ b/src/tests/quadrature/CMakeLists.txt
@@ -1,2 +1,3 @@
 ADDTEST(trapz)
 ADDTEST(simps)
+ADDTEST(gauss)
diff --git a/src/tests/quadrature/Makefile.manual b/src/tests/quadrature/Makefile.manual
index 83061a7bf..024090ea1 100644
--- a/src/tests/quadrature/Makefile.manual
+++ b/src/tests/quadrature/Makefile.manual
@@ -1,3 +1,5 @@
 PROGS_SRC = test_trapz.f90
+PROGS_SRC = test_simps.f90
+PROGS_SRC = test_gauss.f90
 
 include ../Makefile.manual.test.mk
diff --git a/src/tests/quadrature/test_gauss.f90 b/src/tests/quadrature/test_gauss.f90
new file mode 100644
index 000000000..b00e882bb
--- /dev/null
+++ b/src/tests/quadrature/test_gauss.f90
@@ -0,0 +1,34 @@
+program test_gauss_
+    use stdlib_kinds, only: dp
+    use stdlib_error, only: check
+    use stdlib_quadrature , only: gauss_legendre, gauss_legendre_lobatto
+
+    implicit none
+
+    call test_gauss
+    call test_gauss_lobatto
+
+contains
+
+    subroutine test_gauss
+        integer :: i
+        real(dp) :: analytic, numeric
+
+        ! x**2 from -1 to 1
+        analytic = 2.0_dp/3.0_dp
+        do i=2,6
+            block
+                real(dp), dimension(i) :: x,w
+                call gauss_legendre(x,w)
+                numeric = sum(x**2 * w)
+                !print *, i, numeric
+                call check(abs(numeric-analytic) < 2*epsilon(analytic))
+            end block
+        end do
+
+    end subroutine
+
+    subroutine test_gauss_lobatto
+    end subroutine
+
+end program

From 6f06f69f84186054add12a7253de92a1786e0613 Mon Sep 17 00:00:00 2001
From: "Harris M. Snyder" <harris.snyder@gmail.com>
Date: Wed, 23 Jun 2021 17:15:14 -0400
Subject: [PATCH 09/18] test implementation for gaussian quadrature functions

---
 src/tests/quadrature/Makefile.manual |   4 +-
 src/tests/quadrature/test_gauss.f90  | 499 +++++++++++++++++++++++++++
 2 files changed, 501 insertions(+), 2 deletions(-)

diff --git a/src/tests/quadrature/Makefile.manual b/src/tests/quadrature/Makefile.manual
index 024090ea1..f8200913e 100644
--- a/src/tests/quadrature/Makefile.manual
+++ b/src/tests/quadrature/Makefile.manual
@@ -1,5 +1,5 @@
 PROGS_SRC = test_trapz.f90
-PROGS_SRC = test_simps.f90
-PROGS_SRC = test_gauss.f90
+PROGS_SRC += test_simps.f90
+PROGS_SRC += test_gauss.f90
 
 include ../Makefile.manual.test.mk
diff --git a/src/tests/quadrature/test_gauss.f90 b/src/tests/quadrature/test_gauss.f90
index b00e882bb..244550437 100644
--- a/src/tests/quadrature/test_gauss.f90
+++ b/src/tests/quadrature/test_gauss.f90
@@ -14,6 +14,7 @@ subroutine test_gauss
         integer :: i
         real(dp) :: analytic, numeric
 
+        ! test an analytic derivative versus an actual integration
         ! x**2 from -1 to 1
         analytic = 2.0_dp/3.0_dp
         do i=2,6
@@ -26,9 +27,507 @@ subroutine test_gauss
             end block
         end do
 
+        ! test the values of nodes and weights
+        i = 5
+        block 
+            real(dp), dimension(i) :: x,w,xref,wref
+            call gauss_legendre(x,w)
+            xref(1)=-0.90617984593866399_dp
+            xref(2)=-0.53846931010568309_dp
+            xref(3)=0.0_dp
+            xref(4)=0.53846931010568309_dp
+            xref(5)=0.90617984593866399_dp
+
+            wref(1)=0.23692688505618909_dp
+            wref(2)=0.47862867049936647_dp
+            wref(3)=0.56888888888888889_dp
+            wref(4)=0.47862867049936647_dp
+            wref(5)=0.23692688505618909_dp
+
+            call check (all(abs(x-xref) < 2*epsilon(x(1))))
+            call check (all(abs(w-wref) < 2*epsilon(w(1))))
+        end block
+
+        i = 32 
+        block 
+            real(dp), dimension(i) :: x,w,xref,wref
+            call gauss_legendre(x,w)
+
+            xref(1)=-0.99726386184948156_dp
+            xref(2)=-0.98561151154526834_dp
+            xref(3)=-0.96476225558750643_dp
+            xref(4)=-0.93490607593773969_dp
+            xref(5)=-0.89632115576605212_dp
+            xref(6)=-0.84936761373256997_dp
+            xref(7)=-0.79448379596794241_dp
+            xref(8)=-0.73218211874028968_dp
+            xref(9)=-0.66304426693021520_dp
+            xref(10)=-0.58771575724076233_dp
+            xref(11)=-0.50689990893222939_dp
+            xref(12)=-0.42135127613063535_dp
+            xref(13)=-0.33186860228212765_dp
+            xref(14)=-0.23928736225213707_dp
+            xref(15)=-0.14447196158279649_dp
+            xref(16)=-0.048307665687738316_dp
+            xref(17)=0.048307665687738316_dp
+            xref(18)=0.14447196158279649_dp
+            xref(19)=0.23928736225213707_dp
+            xref(20)=0.33186860228212765_dp
+            xref(21)=0.42135127613063535_dp
+            xref(22)=0.50689990893222939_dp
+            xref(23)=0.58771575724076233_dp
+            xref(24)=0.66304426693021520_dp
+            xref(25)=0.73218211874028968_dp
+            xref(26)=0.79448379596794241_dp
+            xref(27)=0.84936761373256997_dp
+            xref(28)=0.89632115576605212_dp
+            xref(29)=0.93490607593773969_dp
+            xref(30)=0.96476225558750643_dp
+            xref(31)=0.98561151154526834_dp
+            xref(32)=0.99726386184948156_dp
+
+
+            wref(1)=0.0070186100094700966_dp
+            wref(2)=0.016274394730905671_dp
+            wref(3)=0.025392065309262059_dp
+            wref(4)=0.034273862913021433_dp
+            wref(5)=0.042835898022226681_dp
+            wref(6)=0.050998059262376176_dp
+            wref(7)=0.058684093478535547_dp
+            wref(8)=0.065822222776361847_dp
+            wref(9)=0.072345794108848506_dp
+            wref(10)=0.078193895787070306_dp
+            wref(11)=0.083311924226946755_dp
+            wref(12)=0.087652093004403811_dp
+            wref(13)=0.091173878695763885_dp
+            wref(14)=0.093844399080804566_dp
+            wref(15)=0.095638720079274859_dp
+            wref(16)=0.096540088514727801_dp
+            wref(17)=0.096540088514727801_dp
+            wref(18)=0.095638720079274859_dp
+            wref(19)=0.093844399080804566_dp
+            wref(20)=0.091173878695763885_dp
+            wref(21)=0.087652093004403811_dp
+            wref(22)=0.083311924226946755_dp
+            wref(23)=0.078193895787070306_dp
+            wref(24)=0.072345794108848506_dp
+            wref(25)=0.065822222776361847_dp
+            wref(26)=0.058684093478535547_dp
+            wref(27)=0.050998059262376176_dp
+            wref(28)=0.042835898022226681_dp
+            wref(29)=0.034273862913021433_dp
+            wref(30)=0.025392065309262059_dp
+            wref(31)=0.016274394730905671_dp
+            wref(32)=0.0070186100094700966_dp
+
+            call check (all(abs(x-xref) < 2*epsilon(x(1))))
+            call check (all(abs(w-wref) < 2*epsilon(w(1))))
+        end block
+
+
+        i = 64
+        block 
+            real(dp), dimension(i) :: x,w,xref,wref
+            call gauss_legendre(x,w)
+
+
+            xref(1)=-0.99930504173577214_dp
+            xref(2)=-0.99634011677195528_dp
+            xref(3)=-0.99101337147674432_dp
+            xref(4)=-0.98333625388462596_dp
+            xref(5)=-0.97332682778991096_dp
+            xref(6)=-0.96100879965205372_dp
+            xref(7)=-0.94641137485840282_dp
+            xref(8)=-0.92956917213193958_dp
+            xref(9)=-0.91052213707850281_dp
+            xref(10)=-0.88931544599511411_dp
+            xref(11)=-0.86599939815409282_dp
+            xref(12)=-0.84062929625258036_dp
+            xref(13)=-0.81326531512279756_dp
+            xref(14)=-0.78397235894334141_dp
+            xref(15)=-0.75281990726053190_dp
+            xref(16)=-0.71988185017161083_dp
+            xref(17)=-0.68523631305423324_dp
+            xref(18)=-0.64896547125465734_dp
+            xref(19)=-0.61115535517239325_dp
+            xref(20)=-0.57189564620263403_dp
+            xref(21)=-0.53127946401989455_dp
+            xref(22)=-0.48940314570705296_dp
+            xref(23)=-0.44636601725346409_dp
+            xref(24)=-0.40227015796399160_dp
+            xref(25)=-0.35722015833766812_dp
+            xref(26)=-0.31132287199021096_dp
+            xref(27)=-0.26468716220876742_dp
+            xref(28)=-0.21742364374000708_dp
+            xref(29)=-0.16964442042399282_dp
+            xref(30)=-0.12146281929612055_dp
+            xref(31)=-0.072993121787799039_dp
+            xref(32)=-0.024350292663424433_dp
+            xref(33)=0.024350292663424433_dp
+            xref(34)=0.072993121787799039_dp
+            xref(35)=0.12146281929612055_dp
+            xref(36)=0.16964442042399282_dp
+            xref(37)=0.21742364374000708_dp
+            xref(38)=0.26468716220876742_dp
+            xref(39)=0.31132287199021096_dp
+            xref(40)=0.35722015833766812_dp
+            xref(41)=0.40227015796399160_dp
+            xref(42)=0.44636601725346409_dp
+            xref(43)=0.48940314570705296_dp
+            xref(44)=0.53127946401989455_dp
+            xref(45)=0.57189564620263403_dp
+            xref(46)=0.61115535517239325_dp
+            xref(47)=0.64896547125465734_dp
+            xref(48)=0.68523631305423324_dp
+            xref(49)=0.71988185017161083_dp
+            xref(50)=0.75281990726053190_dp
+            xref(51)=0.78397235894334141_dp
+            xref(52)=0.81326531512279756_dp
+            xref(53)=0.84062929625258036_dp
+            xref(54)=0.86599939815409282_dp
+            xref(55)=0.88931544599511411_dp
+            xref(56)=0.91052213707850281_dp
+            xref(57)=0.92956917213193958_dp
+            xref(58)=0.94641137485840282_dp
+            xref(59)=0.96100879965205372_dp
+            xref(60)=0.97332682778991096_dp
+            xref(61)=0.98333625388462596_dp
+            xref(62)=0.99101337147674432_dp
+            xref(63)=0.99634011677195528_dp
+            xref(64)=0.99930504173577214_dp
+
+
+            wref(1)=0.0017832807216964329_dp
+            wref(2)=0.0041470332605624676_dp
+            wref(3)=0.0065044579689783629_dp
+            wref(4)=0.0088467598263639477_dp
+            wref(5)=0.011168139460131129_dp
+            wref(6)=0.013463047896718643_dp
+            wref(7)=0.015726030476024719_dp
+            wref(8)=0.017951715775697343_dp
+            wref(9)=0.020134823153530209_dp
+            wref(10)=0.022270173808383254_dp
+            wref(11)=0.024352702568710873_dp
+            wref(12)=0.026377469715054659_dp
+            wref(13)=0.028339672614259483_dp
+            wref(14)=0.030234657072402479_dp
+            wref(15)=0.032057928354851554_dp
+            wref(16)=0.033805161837141609_dp
+            wref(17)=0.035472213256882384_dp
+            wref(18)=0.037055128540240046_dp
+            wref(19)=0.038550153178615629_dp
+            wref(20)=0.039953741132720341_dp
+            wref(21)=0.041262563242623529_dp
+            wref(22)=0.042473515123653589_dp
+            wref(23)=0.043583724529323453_dp
+            wref(24)=0.044590558163756563_dp
+            wref(25)=0.045491627927418144_dp
+            wref(26)=0.046284796581314417_dp
+            wref(27)=0.046968182816210017_dp
+            wref(28)=0.047540165714830309_dp
+            wref(29)=0.047999388596458308_dp
+            wref(30)=0.048344762234802957_dp
+            wref(31)=0.048575467441503427_dp
+            wref(32)=0.048690957009139720_dp
+            wref(33)=0.048690957009139720_dp
+            wref(34)=0.048575467441503427_dp
+            wref(35)=0.048344762234802957_dp
+            wref(36)=0.047999388596458308_dp
+            wref(37)=0.047540165714830309_dp
+            wref(38)=0.046968182816210017_dp
+            wref(39)=0.046284796581314417_dp
+            wref(40)=0.045491627927418144_dp
+            wref(41)=0.044590558163756563_dp
+            wref(42)=0.043583724529323453_dp
+            wref(43)=0.042473515123653589_dp
+            wref(44)=0.041262563242623529_dp
+            wref(45)=0.039953741132720341_dp
+            wref(46)=0.038550153178615629_dp
+            wref(47)=0.037055128540240046_dp
+            wref(48)=0.035472213256882384_dp
+            wref(49)=0.033805161837141609_dp
+            wref(50)=0.032057928354851554_dp
+            wref(51)=0.030234657072402479_dp
+            wref(52)=0.028339672614259483_dp
+            wref(53)=0.026377469715054659_dp
+            wref(54)=0.024352702568710873_dp
+            wref(55)=0.022270173808383254_dp
+            wref(56)=0.020134823153530209_dp
+            wref(57)=0.017951715775697343_dp
+            wref(58)=0.015726030476024719_dp
+            wref(59)=0.013463047896718643_dp
+            wref(60)=0.011168139460131129_dp
+            wref(61)=0.0088467598263639477_dp
+            wref(62)=0.0065044579689783629_dp
+            wref(63)=0.0041470332605624676_dp
+            wref(64)=0.0017832807216964329_dp
+            
+            call check (all(abs(x-xref) < 2*epsilon(x(1))))
+            call check (all(abs(w-wref) < 2*epsilon(w(1))))
+        end block
+
+
+
     end subroutine
 
     subroutine test_gauss_lobatto
+        integer :: i
+        real(dp) :: analytic, numeric
+
+        ! test an analytic derivative versus an actual integration
+        ! x**2 from -1 to 1
+        analytic = 2.0_dp/3.0_dp
+        do i=4,6 ! lobatto quadrature is less accurate for low i, so omit checking at i=2,3
+            block
+                real(dp), dimension(i) :: x,w
+                call gauss_legendre_lobatto(x,w)
+                numeric = sum(x**2 * w)
+                !print *, i, numeric
+                call check(abs(numeric-analytic) < 2*epsilon(analytic))
+            end block
+        end do
+
+
+        ! test the values of nodes and weights
+        i = 5
+        block 
+            real(dp), dimension(i) :: x,w,xref,wref
+            call gauss_legendre_lobatto(x,w)
+
+
+            xref(1)=-1.0000000000000000_dp
+            xref(2)=-0.65465367070797714_dp
+            xref(3)=0.0_dp
+            xref(4)=0.65465367070797714_dp
+            xref(5)=1.0000000000000000_dp
+
+            wref(1)=0.10000000000000000_dp
+            wref(2)=0.54444444444444444_dp
+            wref(3)=0.71111111111111111_dp
+            wref(4)=0.54444444444444444_dp
+            wref(5)=0.10000000000000000_dp
+
+
+            call check (all(abs(x-xref) < 2*epsilon(x(1))))
+            call check (all(abs(w-wref) < 2*epsilon(w(1))))
+        end block
+
+        i = 32 
+        block 
+            real(dp), dimension(i) :: x,w,xref,wref
+            call gauss_legendre_lobatto(x,w)
+
+            xref(1)=-1.0000000000000000_dp
+            xref(2)=-0.99260893397276136_dp
+            xref(3)=-0.97529469048270923_dp
+            xref(4)=-0.94828483841723238_dp
+            xref(5)=-0.91184993906373190_dp
+            xref(6)=-0.86635247601267552_dp
+            xref(7)=-0.81224473177744234_dp
+            xref(8)=-0.75006449393667480_dp
+            xref(9)=-0.68042975561555082_dp
+            xref(10)=-0.60403258714842113_dp
+            xref(11)=-0.52163226288156529_dp
+            xref(12)=-0.43404771720184694_dp
+            xref(13)=-0.34214940653888149_dp
+            xref(14)=-0.24685065885020530_dp
+            xref(15)=-0.14909859681364749_dp
+            xref(16)=-0.049864725046593252_dp
+            xref(17)=0.049864725046593252_dp
+            xref(18)=0.14909859681364749_dp
+            xref(19)=0.24685065885020530_dp
+            xref(20)=0.34214940653888149_dp
+            xref(21)=0.43404771720184694_dp
+            xref(22)=0.52163226288156529_dp
+            xref(23)=0.60403258714842113_dp
+            xref(24)=0.68042975561555082_dp
+            xref(25)=0.75006449393667480_dp
+            xref(26)=0.81224473177744234_dp
+            xref(27)=0.86635247601267552_dp
+            xref(28)=0.91184993906373190_dp
+            xref(29)=0.94828483841723238_dp
+            xref(30)=0.97529469048270923_dp
+            xref(31)=0.99260893397276136_dp
+            xref(32)=1.0000000000000000_dp
+
+            wref(1)=0.0020161290322580645_dp
+            wref(2)=0.012398106501373844_dp
+            wref(3)=0.022199552889291965_dp
+            wref(4)=0.031775135410915466_dp
+            wref(5)=0.041034201586062723_dp
+            wref(6)=0.049885271336221207_dp
+            wref(7)=0.058240497248055870_dp
+            wref(8)=0.066016877257154544_dp
+            wref(9)=0.073137139602679033_dp
+            wref(10)=0.079530525692106252_dp
+            wref(11)=0.085133497949668231_dp
+            wref(12)=0.089890372957357833_dp
+            wref(13)=0.093753875546813814_dp
+            wref(14)=0.096685608948002601_dp
+            wref(15)=0.098656436540761777_dp
+            wref(16)=0.099646771501276778_dp
+            wref(17)=0.099646771501276778_dp
+            wref(18)=0.098656436540761777_dp
+            wref(19)=0.096685608948002601_dp
+            wref(20)=0.093753875546813814_dp
+            wref(21)=0.089890372957357833_dp
+            wref(22)=0.085133497949668231_dp
+            wref(23)=0.079530525692106252_dp
+            wref(24)=0.073137139602679033_dp
+            wref(25)=0.066016877257154544_dp
+            wref(26)=0.058240497248055870_dp
+            wref(27)=0.049885271336221207_dp
+            wref(28)=0.041034201586062723_dp
+            wref(29)=0.031775135410915466_dp
+            wref(30)=0.022199552889291965_dp
+            wref(31)=0.012398106501373844_dp
+            wref(32)=0.0020161290322580645_dp
+            
+            call check (all(abs(x-xref) < 2*epsilon(x(1))))
+            call check (all(abs(w-wref) < 2*epsilon(w(1))))
+        end block
+
+
+        i = 64
+        block 
+            real(dp), dimension(i) :: x,w,xref,wref
+            call gauss_legendre_lobatto(x,w)
+
+
+            xref(1)=-1.0000000000000000_dp
+            xref(2)=-0.99817987150216322_dp
+            xref(3)=-0.99390272670305729_dp
+            xref(4)=-0.98719267660274024_dp
+            xref(5)=-0.97806666283139607_dp
+            xref(6)=-0.96654711036909923_dp
+            xref(7)=-0.95266223578866292_dp
+            xref(8)=-0.93644602747563416_dp
+            xref(9)=-0.91793817351028163_dp
+            xref(10)=-0.89718396784585004_dp
+            xref(11)=-0.87423420065762749_dp
+            xref(12)=-0.84914503454299099_dp
+            xref(13)=-0.82197786730751705_dp
+            xref(14)=-0.79279918182620814_dp
+            xref(15)=-0.76168038340811997_dp
+            xref(16)=-0.72869762508883694_dp
+            xref(17)=-0.69393162129070484_dp
+            xref(18)=-0.65746745031297650_dp
+            xref(19)=-0.61939434613843155_dp
+            xref(20)=-0.57980548006771842_dp
+            xref(21)=-0.53879773271680044_dp
+            xref(22)=-0.49647145693605775_dp
+            xref(23)=-0.45293023223158118_dp
+            xref(24)=-0.40828061128985404_dp
+            xref(25)=-0.36263185922626182_dp
+            xref(26)=-0.31609568619562581_dp
+            xref(27)=-0.26878597401917000_dp
+            xref(28)=-0.22081849749695350_dp
+            xref(29)=-0.17231064108779297_dp
+            xref(30)=-0.12338111165002799_dp
+            xref(31)=-0.074149647946115919_dp
+            xref(32)=-0.024736727621958728_dp
+            xref(33)=0.024736727621958728_dp
+            xref(34)=0.074149647946115919_dp
+            xref(35)=0.12338111165002799_dp
+            xref(36)=0.17231064108779297_dp
+            xref(37)=0.22081849749695350_dp
+            xref(38)=0.26878597401917000_dp
+            xref(39)=0.31609568619562581_dp
+            xref(40)=0.36263185922626182_dp
+            xref(41)=0.40828061128985404_dp
+            xref(42)=0.45293023223158118_dp
+            xref(43)=0.49647145693605775_dp
+            xref(44)=0.53879773271680044_dp
+            xref(45)=0.57980548006771842_dp
+            xref(46)=0.61939434613843155_dp
+            xref(47)=0.65746745031297650_dp
+            xref(48)=0.69393162129070484_dp
+            xref(49)=0.72869762508883694_dp
+            xref(50)=0.76168038340811997_dp
+            xref(51)=0.79279918182620814_dp
+            xref(52)=0.82197786730751705_dp
+            xref(53)=0.84914503454299099_dp
+            xref(54)=0.87423420065762749_dp
+            xref(55)=0.89718396784585004_dp
+            xref(56)=0.91793817351028163_dp
+            xref(57)=0.93644602747563416_dp
+            xref(58)=0.95266223578866292_dp
+            xref(59)=0.96654711036909923_dp
+            xref(60)=0.97806666283139607_dp
+            xref(61)=0.98719267660274024_dp
+            xref(62)=0.99390272670305729_dp
+            xref(63)=0.99817987150216322_dp
+            xref(64)=1.0000000000000000_dp
+
+            wref(1)=0.00049603174603174603_dp
+            wref(2)=0.0030560082449124904_dp
+            wref(3)=0.0054960162038171569_dp
+            wref(4)=0.0079212897900466340_dp
+            wref(5)=0.010327002366815328_dp
+            wref(6)=0.012707399197454735_dp
+            wref(7)=0.015056683987961443_dp
+            wref(8)=0.017369116384542182_dp
+            wref(9)=0.019639040723241718_dp
+            wref(10)=0.021860903511518060_dp
+            wref(11)=0.024029268144023827_dp
+            wref(12)=0.026138828614338438_dp
+            wref(13)=0.028184422665848517_dp
+            wref(14)=0.030161044499089451_dp
+            wref(15)=0.032063857057727025_dp
+            wref(16)=0.033888203884125398_dp
+            wref(17)=0.035629620524489486_dp
+            wref(18)=0.037283845459801173_dp
+            wref(19)=0.038846830537807737_dp
+            wref(20)=0.040314750881560237_dp
+            wref(21)=0.041684014250801952_dp
+            wref(22)=0.042951269833601819_dp
+            wref(23)=0.044113416446892471_dp
+            wref(24)=0.045167610125947702_dp
+            wref(25)=0.046111271084289059_dp
+            wref(26)=0.046942090027028316_dp
+            wref(27)=0.047658033802220637_dp
+            wref(28)=0.048257350376414549_dp
+            wref(29)=0.048738573122233185_dp
+            wref(30)=0.049100524407501308_dp
+            wref(31)=0.049342318477139574_dp
+            wref(32)=0.049463363620776646_dp
+            wref(33)=0.049463363620776646_dp
+            wref(34)=0.049342318477139574_dp
+            wref(35)=0.049100524407501308_dp
+            wref(36)=0.048738573122233185_dp
+            wref(37)=0.048257350376414549_dp
+            wref(38)=0.047658033802220637_dp
+            wref(39)=0.046942090027028316_dp
+            wref(40)=0.046111271084289059_dp
+            wref(41)=0.045167610125947702_dp
+            wref(42)=0.044113416446892471_dp
+            wref(43)=0.042951269833601819_dp
+            wref(44)=0.041684014250801952_dp
+            wref(45)=0.040314750881560237_dp
+            wref(46)=0.038846830537807737_dp
+            wref(47)=0.037283845459801173_dp
+            wref(48)=0.035629620524489486_dp
+            wref(49)=0.033888203884125398_dp
+            wref(50)=0.032063857057727025_dp
+            wref(51)=0.030161044499089451_dp
+            wref(52)=0.028184422665848517_dp
+            wref(53)=0.026138828614338438_dp
+            wref(54)=0.024029268144023827_dp
+            wref(55)=0.021860903511518060_dp
+            wref(56)=0.019639040723241718_dp
+            wref(57)=0.017369116384542182_dp
+            wref(58)=0.015056683987961443_dp
+            wref(59)=0.012707399197454735_dp
+            wref(60)=0.010327002366815328_dp
+            wref(61)=0.0079212897900466340_dp
+            wref(62)=0.0054960162038171569_dp
+            wref(63)=0.0030560082449124904_dp
+            wref(64)=0.00049603174603174603_dp
+            
+            call check (all(abs(x-xref) < 2*epsilon(x(1))))
+            call check (all(abs(w-wref) < 2*epsilon(w(1))))
+        end block
+
     end subroutine
 
 end program

From 274f5423aae74af5d8253f56aab580117b929a5f Mon Sep 17 00:00:00 2001
From: "Harris M. Snyder" <harris.snyder@gmail.com>
Date: Wed, 23 Jun 2021 17:46:11 -0400
Subject: [PATCH 10/18] added spec for gaussian quadrature functions

---
 doc/specs/stdlib_quadrature.md | 84 ++++++++++++++++++++++++++++++++++
 1 file changed, 84 insertions(+)

diff --git a/doc/specs/stdlib_quadrature.md b/doc/specs/stdlib_quadrature.md
index 6f0753367..dcf11e24d 100644
--- a/doc/specs/stdlib_quadrature.md
+++ b/doc/specs/stdlib_quadrature.md
@@ -186,3 +186,87 @@ program demo_simps_weights
 ! 64.0
 end program demo_simps_weights
 ```
+
+## `gauss_legendre` - Gauss-Legendre quadrature (a.k.a. Gaussian quadrature) nodes and weights
+
+### Status
+
+Experimental
+
+### Description
+
+Computes Gauss-Legendre quadrature (also known as simply Gaussian quadrature) nodes and weights,
+ for any `N` (number of nodes).
+Using the nodes `x` and weights `w`, you can compute the integral of some function `f` as follows:
+`integral = sum(f(x) * w)`.
+
+Only double precision is supported - if lower precision is required, you must do the appropriate conversion yourself.
+Accuracy has been validated up to N=64 by comparing computed results to tablulated values known to be accurate to machine precision
+(maximum difference from those values is 2 epsilon).
+
+### Syntax
+
+`subroutine [[stdlib_quadrature(module):gauss_legendre(interface)]](x, w [, interval])`
+
+### Arguments
+
+`x`: Shall be a rank-one array of type `real(real64)`. It is an *output* argument, representing the quadrature nodes.
+
+`w`: Shall be a rank-one array of type `real(real64)`, with the same dimension as `x`. 
+It is an *output* argument, representing the quadrature weights.
+
+`interval`: (Optional) Shall be a two-element array of type `real(real64)`. 
+If present, the nodes and weigts are calculated for integration from `interval(1)` to `interval(2)`.
+If not specified, the default integral is -1 to 1.
+
+### Example
+
+```fortran
+integer, parameter :: N = 6
+real(dp), dimension(N) :: x,w
+call gauss_legendre(x,w)
+integral = sum(x**2 * w)
+```
+
+## `gauss_legendre_lobatto` - Gauss-Legendre-Lobatto quadrature nodes and weights
+
+### Status
+
+Experimental
+
+### Description
+
+Computes Gauss-Legendre-Lobatto quadrature nodes and weights,
+ for any `N` (number of nodes).
+Using the nodes `x` and weights `w`, you can compute the integral of some function `f` as follows:
+`integral = sum(f(x) * w)`.
+
+Only double precision is supported - if lower precision is required, you must do the appropriate conversion yourself.
+Accuracy has been validated up to N=64 by comparing computed results to tablulated values known to be accurate to machine precision
+(maximum difference from those values is 2 epsilon).
+
+### Syntax
+
+`subroutine [[stdlib_quadrature(module):gauss_legendre_lobatto(interface)]](x, w [, interval])`
+
+### Arguments
+
+`x`: Shall be a rank-one array of type `real(real64)`. It is an *output* argument, representing the quadrature nodes.
+
+`w`: Shall be a rank-one array of type `real(real64)`, with the same dimension as `x`. 
+It is an *output* argument, representing the quadrature weights.
+
+`interval`: (Optional) Shall be a two-element array of type `real(real64)`. 
+If present, the nodes and weigts are calculated for integration from `interval(1)` to `interval(2)`.
+If not specified, the default integral is -1 to 1.
+
+### Example
+
+```fortran
+integer, parameter :: N = 6
+real(dp), dimension(N) :: x,w
+call gauss_legendre_lobatto(x,w)
+integral = sum(x**2 * w)
+```
+
+

From f6ac41be6158b7b07685b6de56daba773d96dc37 Mon Sep 17 00:00:00 2001
From: "Harris M. Snyder" <harris.snyder@gmail.com>
Date: Wed, 23 Jun 2021 17:53:41 -0400
Subject: [PATCH 11/18] Added spec for special functions

---
 doc/specs/stdlib_specialfunctions.md | 64 ++++++++++++++++++++++++++++
 1 file changed, 64 insertions(+)
 create mode 100644 doc/specs/stdlib_specialfunctions.md

diff --git a/doc/specs/stdlib_specialfunctions.md b/doc/specs/stdlib_specialfunctions.md
new file mode 100644
index 000000000..d3e512a8e
--- /dev/null
+++ b/doc/specs/stdlib_specialfunctions.md
@@ -0,0 +1,64 @@
+---
+title: specialfunctions
+---
+
+# Special functions
+
+[TOC]
+
+## `legendre` - Calculate Legendre polynomials
+
+### Status
+
+Experimental
+
+### Description
+
+Computes the value of the n-th Legendre polynomial at a specified point.
+Currently only 64 bit floating point is supported.
+
+This is an `elemental` function.
+
+### Syntax
+
+`result = [[stdlib_specialfunctions(module):legendre(interface)]](n, x)`
+
+### Arguments
+
+`n`: Shall be a scalar of type `real(real64)`. 
+
+`x`: Shall be a scalar or array (this function is elemental) of type `real(real64)`. 
+
+### Return value
+
+The function result will be the value of the `n`-th Legendre polynomial, evaluated at `x`.
+
+
+
+## `dlegendre` - Calculate first derivatives of Legendre polynomials
+
+### Status
+
+Experimental
+
+### Description
+
+Computes the value of the first derivative of the n-th Legendre polynomial at a specified point.
+Currently only 64 bit floating point is supported.
+
+This is an `elemental` function.
+
+### Syntax
+
+`result = [[stdlib_specialfunctions(module):dlegendre(interface)]](n, x)`
+
+### Arguments
+
+`n`: Shall be a scalar of type `real(real64)`. 
+
+`x`: Shall be a scalar or array (this function is elemental) of type `real(real64)`. 
+
+### Return value
+
+The function result will be the value of the first derivative of the `n`-th Legendre polynomial, evaluated at `x`.
+

From 4dbd769327efe70c06f0cee17d91cf7337d5e69e Mon Sep 17 00:00:00 2001
From: Harris Snyder <harris.snyder@gmail.com>
Date: Thu, 24 Jun 2021 10:05:33 -0400
Subject: [PATCH 12/18] Update doc/specs/stdlib_quadrature.md

Co-authored-by: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com>
---
 doc/specs/stdlib_quadrature.md | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/doc/specs/stdlib_quadrature.md b/doc/specs/stdlib_quadrature.md
index dcf11e24d..259065200 100644
--- a/doc/specs/stdlib_quadrature.md
+++ b/doc/specs/stdlib_quadrature.md
@@ -206,7 +206,7 @@ Accuracy has been validated up to N=64 by comparing computed results to tablulat
 
 ### Syntax
 
-`subroutine [[stdlib_quadrature(module):gauss_legendre(interface)]](x, w [, interval])`
+`subroutine [[stdlib_quadrature(module):gauss_legendre(interface)]] (x, w[, interval])`
 
 ### Arguments
 
@@ -269,4 +269,3 @@ call gauss_legendre_lobatto(x,w)
 integral = sum(x**2 * w)
 ```
 
-

From 259482a7de8d89b992d28a0121708227a4726c7a Mon Sep 17 00:00:00 2001
From: Harris Snyder <harris.snyder@gmail.com>
Date: Thu, 24 Jun 2021 10:05:46 -0400
Subject: [PATCH 13/18] Update doc/specs/stdlib_quadrature.md

Co-authored-by: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com>
---
 doc/specs/stdlib_quadrature.md | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/doc/specs/stdlib_quadrature.md b/doc/specs/stdlib_quadrature.md
index 259065200..1cbc4d16d 100644
--- a/doc/specs/stdlib_quadrature.md
+++ b/doc/specs/stdlib_quadrature.md
@@ -247,7 +247,7 @@ Accuracy has been validated up to N=64 by comparing computed results to tablulat
 
 ### Syntax
 
-`subroutine [[stdlib_quadrature(module):gauss_legendre_lobatto(interface)]](x, w [, interval])`
+`subroutine [[stdlib_quadrature(module):gauss_legendre_lobatto(interface)]] (x, w[, interval])`
 
 ### Arguments
 
@@ -268,4 +268,3 @@ real(dp), dimension(N) :: x,w
 call gauss_legendre_lobatto(x,w)
 integral = sum(x**2 * w)
 ```
-

From 722598731386e1ed87db024d24aa86d57978bf93 Mon Sep 17 00:00:00 2001
From: Harris Snyder <harris.snyder@gmail.com>
Date: Thu, 24 Jun 2021 10:05:54 -0400
Subject: [PATCH 14/18] Update doc/specs/stdlib_specialfunctions.md

Co-authored-by: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com>
---
 doc/specs/stdlib_specialfunctions.md | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/doc/specs/stdlib_specialfunctions.md b/doc/specs/stdlib_specialfunctions.md
index d3e512a8e..f9962faa2 100644
--- a/doc/specs/stdlib_specialfunctions.md
+++ b/doc/specs/stdlib_specialfunctions.md
@@ -50,7 +50,7 @@ This is an `elemental` function.
 
 ### Syntax
 
-`result = [[stdlib_specialfunctions(module):dlegendre(interface)]](n, x)`
+`result = [[stdlib_specialfunctions(module):dlegendre(interface)]] (n, x)`
 
 ### Arguments
 
@@ -61,4 +61,3 @@ This is an `elemental` function.
 ### Return value
 
 The function result will be the value of the first derivative of the `n`-th Legendre polynomial, evaluated at `x`.
-

From 6647e9b23c17bd090b37e5f154d44a5b0fc3190b Mon Sep 17 00:00:00 2001
From: Harris Snyder <harris.snyder@gmail.com>
Date: Thu, 24 Jun 2021 10:06:01 -0400
Subject: [PATCH 15/18] Update doc/specs/stdlib_specialfunctions.md

Co-authored-by: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com>
---
 doc/specs/stdlib_specialfunctions.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/specs/stdlib_specialfunctions.md b/doc/specs/stdlib_specialfunctions.md
index f9962faa2..15f27577c 100644
--- a/doc/specs/stdlib_specialfunctions.md
+++ b/doc/specs/stdlib_specialfunctions.md
@@ -21,7 +21,7 @@ This is an `elemental` function.
 
 ### Syntax
 
-`result = [[stdlib_specialfunctions(module):legendre(interface)]](n, x)`
+`result = [[stdlib_specialfunctions(module):legendre(interface)]] (n, x)`
 
 ### Arguments
 

From f3d43b358c507a77f7f708d0af41c85695047e14 Mon Sep 17 00:00:00 2001
From: "Harris M. Snyder" <harris.snyder@gmail.com>
Date: Thu, 24 Jun 2021 10:09:53 -0400
Subject: [PATCH 16/18] made spec examples (guassian quadrature) standalone
 programs

---
 doc/specs/stdlib_quadrature.md | 26 ++++++++++++++++++--------
 1 file changed, 18 insertions(+), 8 deletions(-)

diff --git a/doc/specs/stdlib_quadrature.md b/doc/specs/stdlib_quadrature.md
index 1cbc4d16d..6b6be64f3 100644
--- a/doc/specs/stdlib_quadrature.md
+++ b/doc/specs/stdlib_quadrature.md
@@ -222,10 +222,15 @@ If not specified, the default integral is -1 to 1.
 ### Example
 
 ```fortran
-integer, parameter :: N = 6
-real(dp), dimension(N) :: x,w
-call gauss_legendre(x,w)
-integral = sum(x**2 * w)
+program integrate
+	use iso_fortran_env, dp => real64
+	implict none
+
+	integer, parameter :: N = 6
+	real(dp), dimension(N) :: x,w
+	call gauss_legendre(x,w)
+	print *, "integral of x**2 from -1 to 1 is",  sum(x**2 * w)
+end program
 ```
 
 ## `gauss_legendre_lobatto` - Gauss-Legendre-Lobatto quadrature nodes and weights
@@ -263,8 +268,13 @@ If not specified, the default integral is -1 to 1.
 ### Example
 
 ```fortran
-integer, parameter :: N = 6
-real(dp), dimension(N) :: x,w
-call gauss_legendre_lobatto(x,w)
-integral = sum(x**2 * w)
+program integrate
+	use iso_fortran_env, dp => real64
+	implict none
+
+	integer, parameter :: N = 6
+	real(dp), dimension(N) :: x,w
+	call gauss_legendre_lobatto(x,w)
+	print *, "integral of x**2 from -1 to 1 is",  sum(x**2 * w)
+end program
 ```

From 6f9d7e48a36e8a8902fd54ab2b183ddc137ac08d Mon Sep 17 00:00:00 2001
From: Harris Snyder <harris.snyder@gmail.com>
Date: Thu, 24 Jun 2021 10:20:36 -0400
Subject: [PATCH 17/18] Update doc/specs/stdlib_quadrature.md

Co-authored-by: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com>
---
 doc/specs/stdlib_quadrature.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/specs/stdlib_quadrature.md b/doc/specs/stdlib_quadrature.md
index 6b6be64f3..669e42318 100644
--- a/doc/specs/stdlib_quadrature.md
+++ b/doc/specs/stdlib_quadrature.md
@@ -224,7 +224,7 @@ If not specified, the default integral is -1 to 1.
 ```fortran
 program integrate
 	use iso_fortran_env, dp => real64
-	implict none
+	implicit none
 
 	integer, parameter :: N = 6
 	real(dp), dimension(N) :: x,w

From 5c5357972213f526b7f4a118e13db3f43b26aa7b Mon Sep 17 00:00:00 2001
From: Harris Snyder <harris.snyder@gmail.com>
Date: Thu, 24 Jun 2021 10:20:42 -0400
Subject: [PATCH 18/18] Update doc/specs/stdlib_quadrature.md

Co-authored-by: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com>
---
 doc/specs/stdlib_quadrature.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/specs/stdlib_quadrature.md b/doc/specs/stdlib_quadrature.md
index 669e42318..d4d39dfcf 100644
--- a/doc/specs/stdlib_quadrature.md
+++ b/doc/specs/stdlib_quadrature.md
@@ -270,7 +270,7 @@ If not specified, the default integral is -1 to 1.
 ```fortran
 program integrate
 	use iso_fortran_env, dp => real64
-	implict none
+	implicit none
 
 	integer, parameter :: N = 6
 	real(dp), dimension(N) :: x,w