diff --git a/Project.toml b/Project.toml
index 6dd26c8..fba26e4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "LegendrePolynomials"
 uuid = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
diff --git a/docs/Project.toml b/docs/Project.toml
index 7f356b3..5e2e83c 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,10 +1,8 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
 HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97"
 LegendrePolynomials = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
 
 [compat]
 Documenter = "0.26"
-DualNumbers = "0.6"
 HyperDualNumbers = "4"
diff --git a/docs/src/derivatives.md b/docs/src/derivatives.md
index 94cf8e1..dec325a 100644
--- a/docs/src/derivatives.md
+++ b/docs/src/derivatives.md
@@ -1,11 +1,31 @@
 # Derivatives of Legendre Polynomials
 
+## Analytical recursive approach
+
+The Bonnet's recursion formula 
+
+```math
+P_\ell(x) = \left((2\ell-1) x P_{\ell-1}(x) - (\ell-1)P_{\ell - 2}(x)\right)/\ell
+```
+
+may be differentiated an arbitrary number of times analytically to obtain recursion relations for higher derivatives: 
+
+```math
+\frac{d^n P_\ell(x)}{dx^n} = \frac{(2\ell-1)}{\ell} \left(x \frac{d^n P_{\ell-1}(x)}{dx^n} + 
+n \frac{d^{(n-1)} P_{\ell-1}(x)}{dx^{(n-1)}} \right) - \frac{(\ell-1)}{\ell} \frac{d^n P_{\ell-2}(x)}{dx^n}
+```
+
+This provides a simultaneous recursion relation in ``\ell`` as well as ``n``, solving which we may obtain derivatives up to any order. This is the approach used in this package to compute the derivatives of Legendre polynomials.
+
+## Automatic diferentiation
+
 The Julia automatic differentiation framework may be used to compute the derivatives of Legendre polynomials alongside their values. Since the defintions of the polynomials are completely general, they may be called with dual or hyperdual numbers as arguments to evaluate derivarives in one go. 
 We demonstrate one example of this using the package [`HyperDualNumbers.jl`](https://github.com/JuliaDiff/HyperDualNumbers.jl) v4:
 
 ```@meta
 DocTestSetup = quote
 	using LegendrePolynomials
+  using HyperDualNumbers
 end
 ```
 
@@ -84,51 +104,4 @@ Pl_dPl_d2Pl (generic function with 1 method)
 
 julia> Pl_dPl_d2Pl(0.5, lmax = 3)
 ([1.0, 0.5, -0.125, -0.4375], [0.0, 1.0, 1.5, 0.375], [0.0, 0.0, 3.0, 7.5])
-```
-
-# Analytical approach for higher derivatives
-
-Legendre polynomials satisfy the differential equation 
-
-```math
-\frac{d}{dx}\left[(1-x^2)\frac{d P_n}{dx} \right] + n(n+1) P_n(x) = 0
-```
-
-We may rearrange the terms to obtain 
-
-```math
-\frac{d^2 P_n}{dx^2} = \frac{1}{(1-x^2)}\left( 2x \frac{d P_n(x)}{dx} - n(n+1)P_n{x} \right)
-```
-
-We may therefore compute the second derivative from the function and its first derivative. Higher derivatives may further be computed in terms of the lower ones.
-
-We demonstrate the second-derivative computation using the package [`DualNumbers.jl`](https://github.com/JuliaDiff/DualNumbers.jl) v0.5: 
-
-```jldoctest dual
-julia> using DualNumbers
-
-julia> x = 0.5;
-
-julia> xd = Dual(x, one(x));
-
-julia> d2Pl(x, P, dP, n) = (2x * dP - n*(n+1) * P)/(1 - x^2);
-
-julia> function d2Pl(x, n)
-           xd = Dual(x, one(x))
-           y = Pl(xd, n)
-           P, dP = realpart(y), dualpart(y)
-           d2Pl(x, P, dP, n)
-       end;
-
-julia> d2Pl(x, 20)
-32.838787646905985
-```
-
-We may check that this matches the result obtained using `HyperDualNumbers`:
-
-```jldoctest hyperdual
-julia> ε₁ε₂part(Pl(xh, 20))
-32.838787646905985
-```
-
-Unfortunately at this point, higher derivatives need to be evaluated analytically and expresed in terms of lower derivatives.
\ No newline at end of file
+```
\ No newline at end of file
diff --git a/docs/src/index.md b/docs/src/index.md
index bd31d6c..6b7fd75 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -7,7 +7,7 @@ end
 
 # Introduction
 
-Compute [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials) using a 3-term recursion relation (Bonnet’s recursion formula).
+Compute [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials) and their derivatives using a 3-term recursion relation (Bonnet’s recursion formula).
 
 ```math
 P_\ell(x) = \left((2\ell-1) x P_{\ell-1}(x) - (\ell-1)P_{\ell - 2}(x)\right)/\ell
@@ -19,10 +19,12 @@ Currently this package evaluates the standard polynomials that satisfy ``P_\ell(
 \int_{-1}^1 P_m(x) P_n(x) dx = \frac{2}{2n+1} \delta_{mn}.
 ```
 
-There are two main functions: 
+There are four main functions: 
 
 * [`Pl(x,l)`](@ref Pl): this evaluates the Legendre polynomial for a given degree `l` at the argument `x`. The argument needs to satisfy `-1 <= x <= 1`.
 * [`collectPl(x; lmax)`](@ref collectPl): this evaluates all the polynomials for `l` lying in `0:lmax` at the argument `x`. As before the argument needs to lie in the domain of validity. Functionally this is equivalent to `Pl.(x, 0:lmax)`, except `collectPl` evaluates the result in one pass, and is therefore faster. There is also the in-place version [`collectPl!`](@ref) that uses a pre-allocated array.
+* [`dnPl(x, l, n)`](@ref dnPl): this evaluates the ``n``-th derivative of the Legendre polynomial ``P_\ell(x)`` at the argument ``x``. The argument needs to satisfy `-1 <= x <= 1`.
+* [`collectdnPl(x; n, lmax)`](@ref collectdnPl): this evaluates the ``n``-th derivative of all the Legendre polynomials for `l = 0:lmax`. There is also an in-place version [`collectdnPl!`](@ref) that uses a pre-allocated array.
 
 # Quick Start
 
@@ -33,7 +35,14 @@ julia> Pl(0.5, 3)
 -0.4375
 ```
 
-Evaluate all the polynomials for `l` in `0:lmax` as 
+Evaluate the `n`th derivative for one `l` as `dnPl(x, l, n)`:
+
+```jldoctest
+julia> dnPl(0.5, 3, 2)
+7.5
+```
+
+Evaluate all the polynomials for `l` in `0:lmax` as `collectPl(x; lmax)`
 
 ```jldoctest
 julia> collectPl(0.5, lmax = 3)
@@ -44,6 +53,19 @@ julia> collectPl(0.5, lmax = 3)
  -0.4375
 ```
 
+Evaluate all the `n`th derivatives as `collectdnPl(x; lmax, n)`:
+
+```jldoctest
+julia> collectdnPl(0.5, lmax = 5, n = 3)
+6-element OffsetArray(::Array{Float64,1}, 0:5) with eltype Float64 with indices 0:5:
+  0.0
+  0.0
+  0.0
+ 15.0
+ 52.5
+ 65.625
+```
+
 # Increase precision
 
 The precision of the result may be changed by using arbitrary-precision types such as `BigFloat`. For example, using `Float64` arguments we obtain
@@ -69,6 +91,16 @@ julia> setprecision(300) do
 0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333317
 ```
 
+This is particularly important to avoid overflow while computing high-order derivatives. For example:
+
+```jldoctest
+julia> dnPl(0.5, 300, 200) # Float64
+NaN
+
+julia> dnPl(big(1)/2, 300, 200) # BigFloat
+1.738632750542319394663553898425873258768856732308227932150592526951212145232716e+499
+```
+
 # Reference
 
 ```@autodocs
diff --git a/src/LegendrePolynomials.jl b/src/LegendrePolynomials.jl
index 64a699f..137dc73 100644
--- a/src/LegendrePolynomials.jl
+++ b/src/LegendrePolynomials.jl
@@ -5,15 +5,24 @@ using OffsetArrays
 export Pl
 export collectPl
 export collectPl!
+export dnPl
+export collectdnPl
+export collectdnPl!
 
 function checkdomain(x)
 	abs(x) > 1 && throw(DomainError(x,
 		"Legendre Polynomials are defined for arguments lying in -1 ⩽ x ⩽ 1"))
 end
 
+assertnonnegative(l) = (l >= 0 || throw(ArgumentError("l must be >= 0, received " * string(l))))
+
 function checksize(arr, lmax)
 	maximum(axes(arr,1)) >= lmax || throw(ArgumentError("array is not large enough to store all values"))
 end
+function checklength(arr, minlength)
+	length(arr) >= minlength || throw(ArgumentError(
+		"array is not large enough to store all values, require a minimum length of " * string(minlength)))
+end
 
 @inline function Pl_recursion(::Type{T}, ℓ, Plm1, Plm2, x) where {T}
 	# relation is valid from ℓ = 1
@@ -21,13 +30,30 @@ end
 	convert(T, Pl)
 end
 
-polytype(x) = typeof(x*x/1)
+@inline function dPl_recursion(::Type{T}, ℓ, n, P_n_lm1, P_nm1_lm1, P_n_lm2, x) where {T}
+	P_n_l = ((2ℓ - 1) * (x * P_n_lm1 + n * P_nm1_lm1) - (ℓ - 1) * P_n_lm2)/ℓ
+	convert(T, P_n_l)
+end
+
+# special cases
+# case 1 : l == n, in which case there are no (l-1,n) and (l-2,n) terms
+@inline function dPl_recursion(::Type{T}, ℓ, n, ::Nothing, P_nm1_lm1, ::Nothing, x) where {T}
+	P_n_l = (2ℓ-1) * P_nm1_lm1
+	convert(T, P_n_l)
+end
+# case 1 : l == n + 1, in which case there's no (l-2,n) term
+@inline function dPl_recursion(::Type{T}, ℓ, n, P_n_lm1, P_nm1_lm1, ::Nothing, x) where {T}
+	P_n_l = ((2ℓ - 1) * (x * P_n_lm1 + n * P_nm1_lm1))/ℓ
+	convert(T, P_n_l)
+end
+
+polytype(x) = typeof(float(x*x))
 
 """
 	LegendrePolynomialIterator(x, [lmax::Integer])
 
-Return an iterator that generates the values of the Legendre polynomials ``P_l(x)`` for the given `x`.
-If `lmax` is specified then only the values of ``P_l(x)`` from `0` to `lmax` are returned.
+Return an iterator that generates the values of the Legendre polynomials ``P_\\ell(x)`` for the given `x`.
+If `lmax` is specified then only the values of ``P_\\ell(x)`` from `0` to `lmax` are returned.
 
 # Examples
 ```jldoctest
@@ -65,10 +91,14 @@ julia> collect(Iterators.take(Iterators.drop(iter, 100), 5)) # evaluate Pl for l
 struct LegendrePolynomialIterator{T, L <: Union{Integer, Nothing}, V}
 	x :: V
 	lmax :: L
+	function LegendrePolynomialIterator{T,L,V}(x::V, lmax::L) where {T, L <: Union{Integer, Nothing}, V}
+		checkdomain(x)
+		new{T,L,V}(x, lmax)		
+	end
 end
 LegendrePolynomialIterator(x) = LegendrePolynomialIterator{polytype(x), Nothing, typeof(x)}(x, nothing)
 function LegendrePolynomialIterator(x, lmax)
-	lmax >= 0 || throw(ArgumentError("degree must be >= 0"))
+	assertnonnegative(lmax)
 	LegendrePolynomialIterator{polytype(x), typeof(lmax), typeof(x)}(x, lmax)
 end
 
@@ -107,7 +137,7 @@ Base.copy(iter::LegendrePolynomialIterator) = typeof(iter)(iter.x, iter.lmax)
 """
 	Pl(x, l::Integer)
 
-Compute the Legendre Polynomial ``P_l(x)`` for the argument `x` and the degree `l`
+Compute the Legendre Polynomial ``P_\\ell(x)`` for the argument `x` and the degree `l`
 
 # Examples
 ```jldoctest
@@ -119,18 +149,101 @@ julia> Pl(0.5, 20)
 ```
 """
 function Pl(x, l::Integer)
-	# Check if x is within limits
-	checkdomain(x)
-
 	iter = LegendrePolynomialIterator(x)
 	d = Iterators.drop(iter, l) # 0 to l-1
 	first(d)
 end
 
+function doublefactorial(T, n)
+	p = convert(T, one(n))
+	for i in n:-2:1
+		p *= convert(T, i)
+	end
+	convert(T, p)
+end
+
+function _checkvalues(x, l, n)
+	assertnonnegative(l)
+	checkdomain(x)
+	n >= 0 || throw(ArgumentError("order of derivative n must be >= 0"))
+end
+
+"""
+	dnPl(x, l::Integer, n::Integer, [cache::AbstractVector])
+
+Compute the ``n``-th derivative ``d^n P_\\ell(x)/dx^n`` of the Legendre polynomial ``P_\\ell(x)``.
+Optionally a pre-allocated vector `cache` may be provided, which must have a minimum length of `l - n + 1` 
+and may be overwritten during the computation.
+
+The order of the derivative `n` must be non-negative. For `n == 0` this function just returns 
+Legendre polynomials.
+
+# Examples
+
+```jldoctest
+julia> dnPl(0.5, 3, 2) # second derivative of P3(x) at x = 0.5
+7.5
+
+julia> dnPl(0.5, 4, 0) == Pl(0.5, 4) # zero-th order derivative == Pl(x)
+true
+```
+"""
+function dnPl(x, l::Integer, n::Integer, 
+	A = begin
+		_checkvalues(x, l, n)
+		# do not allocate A if the value is trivially zero
+		if l < n
+			return zero(polytype(x))
+		end
+		zeros(polytype(x), l - n + 1)
+	end
+	)
+
+	_checkvalues(x, l, n)
+	checklength(A, l - n + 1)
+
+	cache = OffsetArrays.no_offset_view(A)
+
+	# check if the value is trivially zero in case A is provided in the function call
+	if l < n
+		return zero(eltype(cache))
+	end
+	
+	if n == l # may short-circuit this
+		cache[1] = doublefactorial(eltype(cache), 2l-1)
+	else
+		collectPl!(cache, x, lmax = l - n)
+
+		for ni in 1:n
+			# We denote the terms as P_ni_li
+
+			# li == ni
+			P_nim1_nim1 = cache[1]
+			P_ni_ni = dPl_recursion(eltype(cache), ni, ni, nothing, P_nim1_nim1, nothing, x)
+			cache[1] = P_ni_ni
+
+			# li == ni + 1
+			P_nim1_ni = cache[2]
+			P_ni_nip1 = dPl_recursion(eltype(cache), ni + 1, ni, P_ni_ni, P_nim1_ni, nothing, x)
+			cache[2] = P_ni_nip1
+
+			for li in ni+2:min(l, l - n + ni)
+				P_ni_lim2 = cache[li - ni - 1]
+				P_ni_lim1 = cache[li - ni]
+				P_nim1_lim1 = cache[li - ni + 1]
+				P_ni_li = dPl_recursion(eltype(cache), li, ni, P_ni_lim1, P_nim1_lim1, P_ni_lim2, x)
+				cache[li - ni + 1] = P_ni_li
+			end
+		end
+	end
+
+	return cache[l - n + 1]
+end
+
 """
 	collectPl!(v::AbstractVector, x; [lmax::Integer = length(v) - 1])
 
-Compute the Legendre Polynomials ``P_l(x)`` for the argument `x` and all degrees `l` in `0:lmax`, 
+Compute the Legendre Polynomials ``P_\\ell(x)`` for the argument `x` and all degrees `l` in `0:lmax`, 
 and store them in `v`.
 
 At output `v[firstindex(v) + l] == Pl(x,l)`.
@@ -173,8 +286,8 @@ end
 """
 	collectPl(x; lmax::Integer)
 
-Compute the Legendre Polynomial ``P_l(x)`` for the argument `x` and all degrees `l` in `0:lmax`.
-Return an `OffsetArray` `P` with indices `0:lmax`, where `P[l] == Pl(x,l)`
+Compute the Legendre Polynomial ``P_\\ell(x)`` for the argument `x` and all degrees `l` in `0:lmax`.
+Return `P` with indices `0:lmax`, where `P[l] == Pl(x,l)`
 
 # Examples
 ```jldoctest
@@ -189,4 +302,70 @@ julia> collectPl(0.5, lmax = 4)
 """
 collectPl(x; lmax::Integer) = collect(LegendrePolynomialIterator(x, lmax))
 
+"""
+	collectdnPl(x; n::Integer, lmax::Integer)
+
+Compute the ``n``-th derivative of a Legendre Polynomial ``P_\\ell(x)`` for the argument `x` and all degrees `l = 0:lmax`. 
+
+The order of the derivative `n` must be greater than or equal to zero.
+
+Returns `v` with indices `0:lmax`, where `v[l] == dnPl(x, l, n)`.
+
+# Examples
+
+```jldoctest
+julia> collectdnPl(0.5, lmax = 3, n = 2)
+4-element OffsetArray(::Array{Float64,1}, 0:3) with eltype Float64 with indices 0:3:
+ 0.0
+ 0.0
+ 3.0
+ 7.5
+```
+"""
+function collectdnPl(x; lmax::Integer, n::Integer)
+	assertnonnegative(lmax)
+	n >= 0 || throw(ArgumentError("order of derivative n must be >= 0"))
+	v = zeros(polytype(x), 0:lmax)
+	if lmax >= n
+		collectdnPl!(parent(v), x; lmax = lmax, n = n)
+	end
+	v
+end
+
+"""
+	collectdnPl!(v::AbstractVector, x; lmax::Integer, n::Integer)
+
+Compute the ``n``-th derivative of a Legendre Polynomial ``P_\\ell(x)`` for the argument `x` and all degrees `l = 0:lmax`, 
+and store the result in `v`. 
+
+The order of the derivative `n` must be greater than or equal to zero.
+
+At output, `v[l + firstindex(v)] == dnPl(x, l, n)` for `l = 0:lmax`.
+
+# Examples
+
+```jldoctest
+julia> v = zeros(4);
+
+julia> collectdnPl!(v, 0.5, lmax = 3, n = 2)
+4-element Array{Float64,1}:
+ 0.0
+ 0.0
+ 3.0
+ 7.5
+```
+"""
+function collectdnPl!(v, x; lmax::Integer, n::Integer)
+	assertnonnegative(lmax)
+	n >= 0 || throw(ArgumentError("order of derivative n must be >= 0"))
+	checklength(v, lmax + 1)
+	
+	# trivially zero for l < n
+	fill!((@view v[(0:n-1) .+ firstindex(v)]), zero(eltype(v)))
+	# populate the other elements
+	dnPl(x, lmax, n, @view v[(n:lmax) .+ firstindex(v)])
+	
+	v
+end
+
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 3b5f935..5f82cf9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -17,29 +17,38 @@ tohyper(x) = Hyper(x, one(x), one(x), zero(x))
     iter = LegendrePolynomialIterator(0.5, lmax);
     @test length(iter) == lmax + 1
     @test eltype(iter) == Float64
+    @test eltype(typeof(iter)) == Float64
     @test size(iter) == (lmax+1,)
     @test axes(iter) == (0:lmax,)
     @test keys(iter) == 0:lmax
 
     iter2 = copy(iter)
     @test typeof(iter2) == typeof(iter)
-    @test iter2.x == iter.x 
+    @test iter2.x == iter.x
     @test iter2.lmax == iter.lmax
+
+    iter = LegendrePolynomialIterator(0.5);
+    @test Base.IteratorSize(typeof(iter)) == Base.IsInfinite()
+
+    @test_throws DomainError LegendrePolynomialIterator(-1.1, lmax)
+    @test_throws DomainError LegendrePolynomialIterator(-1.1)
 end
 
 @testset "Pl" begin
 	x = 2rand() - 1 
-    P = collectPl(x, lmax = 5)
-    @test P[0] ≈ 1
-    @test P[1] ≈ x
-    @test P[2] ≈ (3x^2 - 1)/2
-    @test P[3] ≈ (5x^3 - 3x)/2
-    @test P[4] ≈ (35x^4 - 30x^2 + 3)/8
-    @test P[5] ≈ (63x^5 - 70x^3 + 15x)/8
+    lmax = 5
+    P = collectPl(x, lmax = lmax)
+    P2 = collectdnPl(x, lmax = lmax, n = 0)
+    @test P[0] == P2[0] == 1
+    @test P[1] == P2[1] == x
+    @test P[2] ≈ P2[2] ≈ (3x^2 - 1)/2
+    @test P[3] ≈ P2[3] ≈ (5x^3 - 3x)/2
+    @test P[4] ≈ P2[4] ≈ (35x^4 - 30x^2 + 3)/8
+    @test P[5] ≈ P2[5] ≈ (63x^5 - 70x^3 + 15x)/8
 
     @testset "x = 0" begin
         for l = 1:2:101
-            @test Pl(0, l) == 0
+            @test Pl(0, l) == dnPl(0, l, 0) == 0
         end
     end
 
@@ -48,12 +57,32 @@ end
 
         P = collectPl(1, lmax = lmax)
         @test all(==(1), P)
+        P2 = collectdnPl(1, lmax = lmax, n = 0)
+        @test all(==(1), P2)
 
         P = collectPl(-1, lmax = lmax)
+        P2 = collectdnPl(-1, lmax = lmax, n = 0)
         for l in axes(P, 1)
-            @test P[l] == (-1)^l
+            @test P[l] == P2[l] == (-1)^l
         end
     end
+
+    @test_throws DomainError collectPl(-2, lmax = lmax)
+    @test_throws ArgumentError collectPl(0, lmax = -1)
+end
+
+@testset "dnPl" begin
+    @test dnPl(0.5, 3, 10) == 0
+    @test dnPl(0.5, 3, 10, zeros(100)) == 0
+    @testset "double factorial" begin
+        @test dnPl(0.5, 2, 2) == 3
+        @test dnPl(0.5, 3, 3) == 15
+        @test dnPl(0.5, 3, 3, zeros(1)) == 15
+    end
+    @test_throws ArgumentError dnPl(0.5, 10, 3, zeros(1))
+    @test_throws ArgumentError dnPl(0.5, -1, 3)
+    @test_throws ArgumentError dnPl(0.5, 1, -3)
+    @test_throws DomainError dnPl(2.5, 1, 3)
 end
 
 @testset "collectPl!" begin
@@ -71,32 +100,56 @@ end
     for l in 0:lmax
         @test v[firstindex(v) + l] == Pl(x, l)
     end
+
+    @test_throws DomainError collectPl!(v, -2, lmax = lmax)
+    @test_throws ArgumentError collectPl!(v, 0, lmax = -1)
+end
+
+@testset "collectdnPl!" begin
+    v = zeros(10);
+    @test_throws DomainError collectdnPl!(v, -2, lmax = 3, n = 2)
+    @test_throws ArgumentError collectdnPl!(v, 0, lmax = -1, n = 3)
+    @test_throws ArgumentError collectdnPl!(v, 0.5, lmax = 1, n = -3)
+    @test_throws ArgumentError collectdnPl!(zeros(1), 0, lmax = 3, n = 3)
+end
+
+@testset "collectdnPl" begin
+    @test_throws DomainError collectdnPl(-2, lmax = 3, n = 2)
+    @test_throws ArgumentError collectdnPl(0, lmax = -1, n = 3)
+    @test_throws ArgumentError collectdnPl(0.5, lmax = 1, n = -3)
 end
 
 @testset "dPl/dx" begin
     x = 2rand() - 1 
     xh = tohyper(x)
-    P = eps1.(collectPl(xh, lmax=5))
-    @test P[0] ≈ 0
-    @test P[1] ≈ 1
-    @test P[2] ≈ 3x
-    @test P[3] ≈ (-3 + 15x^2)/2
-    @test P[4] ≈ 1/8*(-60x + 140x^3)
-    @test P[5] ≈ 1/8*(15 - 210x^2 + 315x^4)
+    lmax = 5
+    P = eps1.(collectPl(xh, lmax = lmax))
+    P2 = collectdnPl(x, lmax = lmax, n = 1)
+    @test P[0] == P2[0] == 0
+    @test P[1] == P2[1] == 1
+    @test P[2] ≈ P2[2] ≈ 3x
+    @test P[3] ≈ P2[3] ≈ (-3 + 15x^2)/2
+    @test P[4] ≈ P2[4] ≈ 1/8*(-60x + 140x^3)
+    @test P[5] ≈ P2[5] ≈ 1/8*(15 - 210x^2 + 315x^4)
 
     @testset "x = ±1" begin
         lmax = 10
-        
-        P = collectPl(tohyper(1), lmax = lmax)
-        dP1 = eps1.(P)
+        x = 1
+        P = collectPl(tohyper(x), lmax = lmax)
+        dP1h = eps1.(P)
+        dP1 = collectdnPl(x, lmax = lmax, n = 1)
+
         for l in axes(dP1,1)
-            @test dP1[l] == l*(l+1)/2
+            @test dP1h[l] == dP1[l] == l*(l+1)/2
         end
 
-        P = collectPl(tohyper(-1), lmax = lmax)
-        dPm1 = eps1.(P)
+        x = -1
+        P = collectPl(tohyper(x), lmax = lmax)
+        dPm1h = eps1.(P)
+        dPm1 = collectdnPl(x, lmax = lmax, n = 1)
+
         for l in axes(dPm1,1)
-            @test dPm1[l] == (-1)^(l+1) * l*(l+1)/2
+            @test dPm1h[l] == dPm1[l] == (-1)^(l+1) * l*(l+1)/2
         end
     end
 end
@@ -104,18 +157,23 @@ end
 @testset "d^2Pl/dx^2" begin
     x = 2rand() - 1 
     xh = tohyper(x)
-    P = eps1eps2.(collectPl(xh, lmax=5))
-    @test P[0] ≈ 0
-    @test P[1] ≈ 0
-    @test P[2] ≈ 3
-    @test P[3] ≈ 15x
-    @test P[4] ≈ 1/8*(-60 + 420x^2)
-    @test P[5] ≈ 1/8*(-420x + 1260x^3)
+    lmax = 5
+    P = eps1eps2.(collectPl(xh, lmax = lmax))
+    P2 = collectdnPl(x, lmax = lmax, n = 2)
+
+    @test P[0] == P2[0] == 0
+    @test P[1] == P2[1] == 0
+
+    @test P[2] ≈ P2[2] ≈ 3
+    @test P[3] ≈ P2[3] ≈ 15x
+    @test P[4] ≈ P2[4] ≈ 1/8*(-60 + 420x^2)
+    @test P[5] ≈ P2[5] ≈ 1/8*(-420x + 1260x^3)
 end
 
 @testset "Askey-Gasper" begin
     x = 0.5
     for lmax=1:100:1_000
         @test sum(collectPl(x, lmax = lmax)) >= 0
+        @test sum(collectdnPl(x, lmax = lmax, n = 0)) >= 0
     end
 end
\ No newline at end of file