diff --git a/README.md b/README.md
index 1d093bc..648ecc7 100644
--- a/README.md
+++ b/README.md
@@ -19,18 +19,18 @@ vx=range(0,1,length=100)
 vy=range(0,2,length=200)
 vz=range(0,3,length=300)
 M=[x^2+y^2+z^2 for x=vx,y=vy,z=vz]
-@show trapz(vx,trapz(vy,trapz(vz,M)));
+@show trapz((vx,vy,vz), M);
 
 ```
 
-    trapz(vx, trapz(vy, trapz(vz, M))) = 28.00030370797026
+    trapz((vx, vy, vz), M) = 28.000303707970264
 
 
 # Benchmarks
 
 
 ```julia
-@benchmark trapz($vx,trapz($vy,trapz($vz,$M)))
+@benchmark trapz(($vx,$vy,$vz),$M)
 ```
 
 
@@ -52,7 +52,7 @@ M=[x^2+y^2+z^2 for x=vx,y=vy,z=vz]
 
 
 ```julia
-@benchmark trapz($vy,trapz($vz,$M))
+@benchmark trapz(($vy, $vz),$M)
 ```
 
 
@@ -100,7 +100,7 @@ This code is optimized in order to perform the integral the fastest over the las
 
 
 ```julia
-@benchmark trapz($vz,trapz($vy,trapz($vx,$M,1),1),1)
+@benchmark trapz(($vz,$vy,$vx),$M,(3,2,1))
 ```
 
 
diff --git a/src/Trapz.jl b/src/Trapz.jl
index 021daa2..7c6b3e1 100644
--- a/src/Trapz.jl
+++ b/src/Trapz.jl
@@ -46,6 +46,25 @@ module Trapz
         trapz(x,PermutedDimsArray(y,bringlast(ntuple(identity,Val(N)),axis)))
     end
 
+    function gen_trapz_impl(xs::Type{T}, M) where T <: NTuple{N} where N
+        ex = :(trapz(xs[$N], M))
+        for i=(N-1):-1:1
+            ex = :(trapz(xs[$i],$ex))
+        end
+        return ex
+    end
+
+    @generated function trapz(xs::NTuple{N}, M) where N
+        return gen_trapz_impl(xs, M)
+    end
+
+    function trapz(xs::NTuple{N}, M::AbstractArray{T,S}, axes::NTuple{N}) where {N, T, S}
+        if S > N
+            axes = ((i for i=1:S if !in(i,axes))..., axes...)
+        end
+        trapz(xs, PermutedDimsArray(M, axes))
+    end
+
     """
         trapz(x,y,axis=End)
         Calculates ∫y[..., i (axis) ,...] dx[i]
@@ -63,7 +82,16 @@ module Trapz
         vx=0:0.01:1
         vy=0:0.01:2
         z=[x^2+y^2 for x=vx, y=vy]
-        trapz(trapz(vy,z),vx)
+        trapz((vx,vy),z) # equivalent to trapz(vx, trapz(vy, z))
+    ```
+            Result ≈ 4/3
+
+        2-D Example in reverse integration order:
+    ```julia
+        vx=0:0.01:1
+        vy=0:0.01:2
+        z=[x^2+y^2 for x=vx, y=vy]
+        trapz((vy,vx),z,(2,1)) # equivalent to trapz(vy, trapz(vx, z, 1), 1)
     ```
             Result ≈ 4/3
     """