diff --git a/julia/MOLE.jl/examples/elliptic1D.jl b/julia/MOLE.jl/examples/elliptic1D.jl
index 0d9583de..c98edb15 100644
--- a/julia/MOLE.jl/examples/elliptic1D.jl
+++ b/julia/MOLE.jl/examples/elliptic1D.jl
@@ -3,7 +3,7 @@ A Julia script that solves the 1D Poisson's equation with Robin boundary conditi
 """
 
 using Plots
-using MOLE
+import MOLE: Operators, BCs
 
 # Domain limits
 west = 0.0
@@ -13,12 +13,12 @@ k = 6 # operator order of accuracy
 m = 2*k + 1 # minimum number of cells to attain desired accuracy
 dx = (east-west)/m # step length
 
-L = lap(k,m,dx)
+L = Operators.lap(k,m,dx)
 
 # Impose Robin boundary condition on laplacian operator
 a = 1.0
 b = 1.0
-L = L + robinBC(k,m,dx,a,b)
+L = L + BCs.robinBC(k,m,dx,a,b)
 
 # 1D staggered grid
 grid = [west; (west+(dx/2)):dx:(east-(dx/2)); east]
diff --git a/julia/MOLE.jl/src/BCs/BCs.jl b/julia/MOLE.jl/src/BCs/BCs.jl
new file mode 100644
index 00000000..0937a8ee
--- /dev/null
+++ b/julia/MOLE.jl/src/BCs/BCs.jl
@@ -0,0 +1,5 @@
+module BCs
+
+include("robinBC.jl")
+
+end # module BCs
\ No newline at end of file
diff --git a/julia/MOLE.jl/src/robinBC.jl b/julia/MOLE.jl/src/BCs/robinBC.jl
similarity index 90%
rename from julia/MOLE.jl/src/robinBC.jl
rename to julia/MOLE.jl/src/BCs/robinBC.jl
index 960ee07f..0fbd5e57 100644
--- a/julia/MOLE.jl/src/robinBC.jl
+++ b/julia/MOLE.jl/src/BCs/robinBC.jl
@@ -4,6 +4,8 @@
     See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details.
 =#
 
+import ..Operators: Operators, grad
+
 """
     robinBC(k, m, dx, a, b)
 
@@ -25,7 +27,7 @@ function robinBC(k::Int, m::Int, dx, a, b)
     B[1,1] = -b
     B[m+2,m+1] = b
 
-    G = grad(k,m,dx)
+    G = Operators.grad(k,m,dx)
 
     BC = A + B*G;
 end
\ No newline at end of file
diff --git a/julia/MOLE.jl/src/MOLE.jl b/julia/MOLE.jl/src/MOLE.jl
index 2c20ae65..5360c41c 100644
--- a/julia/MOLE.jl/src/MOLE.jl
+++ b/julia/MOLE.jl/src/MOLE.jl
@@ -6,11 +6,7 @@
 
 module MOLE
 
-include("divergence.jl")
-include("gradient.jl")
-include("laplacian.jl")
-include("robinBC.jl")
-
-export div, grad, lap, robinBC
+include("Operators/Operators.jl")
+include("BCs/BCs.jl")
 
 end # module MOLE
\ No newline at end of file
diff --git a/julia/MOLE.jl/src/Operators/Operators.jl b/julia/MOLE.jl/src/Operators/Operators.jl
new file mode 100644
index 00000000..20491d46
--- /dev/null
+++ b/julia/MOLE.jl/src/Operators/Operators.jl
@@ -0,0 +1,7 @@
+module Operators
+
+include("gradient.jl")
+include("divergence.jl")
+include("laplacian.jl")
+
+end # module Operators
\ No newline at end of file
diff --git a/julia/MOLE.jl/src/divergence.jl b/julia/MOLE.jl/src/Operators/divergence.jl
similarity index 100%
rename from julia/MOLE.jl/src/divergence.jl
rename to julia/MOLE.jl/src/Operators/divergence.jl
diff --git a/julia/MOLE.jl/src/gradient.jl b/julia/MOLE.jl/src/Operators/gradient.jl
similarity index 100%
rename from julia/MOLE.jl/src/gradient.jl
rename to julia/MOLE.jl/src/Operators/gradient.jl
diff --git a/julia/MOLE.jl/src/laplacian.jl b/julia/MOLE.jl/src/Operators/laplacian.jl
similarity index 100%
rename from julia/MOLE.jl/src/laplacian.jl
rename to julia/MOLE.jl/src/Operators/laplacian.jl
diff --git a/julia/MOLE.jl/test/testDivergence.jl b/julia/MOLE.jl/test/Operators/divergence.jl
similarity index 74%
rename from julia/MOLE.jl/test/testDivergence.jl
rename to julia/MOLE.jl/test/Operators/divergence.jl
index 9a0dd6d6..45c294cc 100644
--- a/julia/MOLE.jl/test/testDivergence.jl
+++ b/julia/MOLE.jl/test/Operators/divergence.jl
@@ -1,8 +1,7 @@
-tol = 1e-10
-
 @testset "Testing divergence for order k=$k" for k=2:2:8
+    tol = 1e-10
     m = 2*k+1
-    D = MOLE.div(k, m, 1/m)
+    D = Operators.div(k, m, 1/m)
     field = ones(m+1,1)
     sol = D*field
     @test norm(sol) < tol
diff --git a/julia/MOLE.jl/test/testGradient.jl b/julia/MOLE.jl/test/Operators/gradient.jl
similarity index 82%
rename from julia/MOLE.jl/test/testGradient.jl
rename to julia/MOLE.jl/test/Operators/gradient.jl
index 3cc0ffd4..d548f1c7 100644
--- a/julia/MOLE.jl/test/testGradient.jl
+++ b/julia/MOLE.jl/test/Operators/gradient.jl
@@ -2,7 +2,7 @@ tol = 1e-10
 
 @testset "Testing gradient for order k=$k" for k=2:2:8
     m = 2*k+1
-    G = MOLE.grad(k, m, 1/m)
+    G = Operators.grad(k, m, 1/m)
     field = ones(m+2, 1)
     sol = G*field
     @test norm(sol) < tol
diff --git a/julia/MOLE.jl/test/testLaplacian.jl b/julia/MOLE.jl/test/Operators/laplacian.jl
similarity index 82%
rename from julia/MOLE.jl/test/testLaplacian.jl
rename to julia/MOLE.jl/test/Operators/laplacian.jl
index 8103589c..7091192d 100644
--- a/julia/MOLE.jl/test/testLaplacian.jl
+++ b/julia/MOLE.jl/test/Operators/laplacian.jl
@@ -2,7 +2,7 @@ tol = 1e-10
 
 @testset "Testing laplacian for order k=$k" for k=2:2:8
     m = 2*k+1
-    L = MOLE.lap(k, m, 1/m)
+    L = Operators.lap(k, m, 1/m)
     field = ones(m+2, 1)
     sol = L*field
     @test norm(sol) < tol
diff --git a/julia/MOLE.jl/test/runtests.jl b/julia/MOLE.jl/test/runtests.jl
index 86b8000a..b2f5274f 100644
--- a/julia/MOLE.jl/test/runtests.jl
+++ b/julia/MOLE.jl/test/runtests.jl
@@ -1,26 +1,17 @@
 using MOLE
 using Test, LinearAlgebra
-
-#tol = 1e-10
+import MOLE: Operators
 
 @testset "Testing MOLE operators" begin
     @testset "Testing 1D divergence" begin
-        include("testDivergence.jl")
+        include("Operators/divergence.jl")
     end
 
     @testset "Testing 1D gradient" begin
-        include("testGradient.jl")
+        include("Operators/gradient.jl")
     end
 
     @testset "Testing 1D laplacian" begin
-        include("testLaplacian.jl")
+        include("Operators/laplacian.jl")
     end
-end
-
-# @testset "Testing Divergence for order k=$k" for k=2:2:8
-#     m = 2*k+1
-#     D = div(k, m, 1/m)
-#     field = ones(m+1,1)
-#     sol = D*field
-#     @test norm(sol) < tol
-# end
\ No newline at end of file
+end
\ No newline at end of file