diff --git a/src/sobol.jl b/src/sobol.jl
index b5e509e63c5dd1ceadeaedab9d3948832302652b..64b455dfd7a8429f9310652cbf14b0f65b07328e 100644
--- a/src/sobol.jl
+++ b/src/sobol.jl
@@ -32,7 +32,16 @@ function sobol(dim::Int64, N::Int64)
         X[j,:], nextseed, MeM  = getNewSobolVector(dim, nextseed, MeM)
     end
 
-    return X
+    return X, nextseed, MeM
+end
+
+function sobol(X, nextseed, MeM)
+    (N,dim) = size(X)
+    for j = 1:N
+        X[j,:], nextseed, MeM  = getNewSobolVector(dim, nextseed, MeM)
+    end
+
+    return X, nextseed, MeM
 end
 
 ## Functions below are inner functions, I have not determined if this is needed, the fastest method or anything else
@@ -106,7 +115,7 @@ end
 function InitSobol(dim::Int64)
     MeM = Mem(dim)
     MeM.seed_save = -1
-#     MeM.v = zeros(40,30) # done in constructor instead
+    #     MeM.v = zeros(40,30) # done in constructor instead
 
     MeM.v[1:40,1] = 1
     MeM.v[3:40,2] = [1, 3, 1, 3, 1, 3, 3, 1,
@@ -188,4 +197,22 @@ function InitSobol(dim::Int64)
     return MeM
 end
 
+"""
+`using Winston`
+Plots the first 4 outputs of sobel(2,512)
+"""
+function test_sobol()
+    X, nextseed, MeM = sobol(2,512)
+    figure(); pp = plot(X[:,1],X[:,2],".b"); hold(true)
+
+    X, nextseed, MeM = sobol(X, nextseed, MeM )
+    plot(X[:,1],X[:,2],".g")
+
+    X, nextseed, MeM = sobol(X, nextseed, MeM )
+    plot(X[:,1],X[:,2],".r")
+
+    X, nextseed, MeM = sobol(X, nextseed, MeM )
+    plot(X[:,1],X[:,2],".m")
 
+    display(pp)
+end
\ No newline at end of file