diff --git a/julia/MOLE.jl/Project.toml b/julia/MOLE.jl/Project.toml
new file mode 100644
index 00000000..e867e80b
--- /dev/null
+++ b/julia/MOLE.jl/Project.toml
@@ -0,0 +1,7 @@
+name = "MOLE"
+uuid = "32da6bfa-ea95-4826-ae25-7537d33b039c"
+version = "0.1.0"
+
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/julia/MOLE.jl/docs/Project.toml b/julia/MOLE.jl/docs/Project.toml
new file mode 100644
index 00000000..dfa65cd1
--- /dev/null
+++ b/julia/MOLE.jl/docs/Project.toml
@@ -0,0 +1,2 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
diff --git a/julia/MOLE.jl/docs/make.jl b/julia/MOLE.jl/docs/make.jl
new file mode 100644
index 00000000..5c8325a2
--- /dev/null
+++ b/julia/MOLE.jl/docs/make.jl
@@ -0,0 +1,9 @@
+push!(LOAD_PATH,"../src/")
+using Documenter, MOLE
+
+makedocs(
+    sitename="MOLE.jl Docs",
+    pages = [
+        "Home" => "index.md"
+    ]
+)
\ No newline at end of file
diff --git a/julia/MOLE.jl/docs/src/index.md b/julia/MOLE.jl/docs/src/index.md
new file mode 100644
index 00000000..cd28af2c
--- /dev/null
+++ b/julia/MOLE.jl/docs/src/index.md
@@ -0,0 +1,58 @@
+# MOLE: Mimetic Operators Library Enhanced
+
+MOLE is a library that implements high-order mimetic operators to solve partial differential equations.
+This site provides documentation for the Julia implementation of MOLE.
+More information can be found on the [main documentation site](https://mole-docs.readthedocs.io/en/main/) and the [GitHub repository](https://github.com/csrc-sdsu/mole).
+
+## Getting Started
+
+MOLE.jl is not yet available in the Julia's package manager. For now, this repository needs to be cloned locally in order to use the library.
+
+## Using MOLE.jl
+
+In order to use the MOLE.jl library, first navigate to the location where the repository has been cloned to. Then, go to the ```mole/julia/MOLE.jl``` sub-directory. From here, you can access the library via the REPL or the command line.
+
+### From the REPL
+
+In ```mole/julia/MOLE.jl```, start Julia using the command ```julia --project=.```. Next, the package lubrary needs to be instantiated and pre-compiled. Enter the ```pkg``` mode by pressing ```]```, then type the commands ```instantiate``` and ```precompile``` (one at a time). This should activate the MOLE.jl package and install the necessary dependencies.
+
+Then, to run a script such as myScript.jl, you can use the following command in the REPL:
+> ```include("path/to/myScript.jl")```
+
+### From the command line
+
+In ```mole/julia/MOLE.jl```, use the following commands to instatiate and precompile the MOLE package:
+
+> ```julia --project=. -e 'using Pkg; Pkg.instantiate()'```
+
+> ```julia --project=. -e 'using Pkg; Pkg.precompile()'```
+
+Then, to run a script such as myScript.jl, use the command
+
+> ```julia --project=. path/to/myScript.jl```
+
+## Functions
+
+### Operators
+
+```@docs
+MOLE.div(k::Int,m::Int,dx)
+MOLE.grad(k::Int,m::Int,dx)
+MOLE.lap(k::Int,m::Int,dx)
+```
+
+### Utilities
+
+```@docs
+MOLE.robinBC(k::Int, m::Int, dx, a, b)
+```
+
+## Examples
+
+The MOLE library contains examples demonstrating how to use the operators, in a broad range of partial differential equations (PDEs). More information on the mathematical content can be found in the [main MOLE documentation](https://mole-docs.readthedocs.io/en/main/examples/index.html).
+
+Currently, the following examples are available in the MOLE Julia package.
+
+- Elliptic Problems
+    - 1D Examples
+        - ```elliptic1D```: A script that solves the 1D Poisson's equation with Robin boundary conditions using mimetic operators.
\ No newline at end of file
diff --git a/julia/MOLE.jl/examples/elliptic1D.jl b/julia/MOLE.jl/examples/elliptic1D.jl
new file mode 100644
index 00000000..0d9583de
--- /dev/null
+++ b/julia/MOLE.jl/examples/elliptic1D.jl
@@ -0,0 +1,39 @@
+"""
+A Julia script that solves the 1D Poisson's equation with Robin boundary conditions using mimetic operators.
+"""
+
+using Plots
+using MOLE
+
+# Domain limits
+west = 0.0
+east = 1.0
+
+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)
+
+# Impose Robin boundary condition on laplacian operator
+a = 1.0
+b = 1.0
+L = L + robinBC(k,m,dx,a,b)
+
+# 1D staggered grid
+grid = [west; (west+(dx/2)):dx:(east-(dx/2)); east]
+
+# RHS
+U = exp.(grid)
+U[1] = 0
+U[end] = 2*exp(1)
+
+U = L\U
+
+# Plot result
+p = scatter(grid, U, label="Approximated")
+plot!(p, grid, exp.(grid), label="Analytical")
+plot!(p, xlabel="x", ylabel="u(x)", title="Poisson's equation with Robin BC")
+display(p)
+println("Press Enter to close the plot and exit.")
+readline()
\ No newline at end of file
diff --git a/julia/MOLE.jl/src/MOLE.jl b/julia/MOLE.jl/src/MOLE.jl
new file mode 100644
index 00000000..a3fcb984
--- /dev/null
+++ b/julia/MOLE.jl/src/MOLE.jl
@@ -0,0 +1,10 @@
+module MOLE
+
+include("divergence.jl")
+include("gradient.jl")
+include("laplacian.jl")
+include("robinBC.jl")
+
+export div, grad, lap, robinBC
+
+end # module MOLE
\ No newline at end of file
diff --git a/julia/MOLE.jl/src/divergence.jl b/julia/MOLE.jl/src/divergence.jl
new file mode 100644
index 00000000..0a3d54cd
--- /dev/null
+++ b/julia/MOLE.jl/src/divergence.jl
@@ -0,0 +1,56 @@
+"""
+    div(k, m, dx)
+
+Returns a m+2 by m+1 one-dimensional mimetic divergence operator.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells
+- `dx`: Step size
+"""
+function div(k::Int,m::Int,dx)
+    if k < 2
+        throw(DomainError(k, "k must be >= 2"))
+    end
+
+    if k % 2 != 0
+        throw(DomainError(k, "k must be an positive even integer"))
+    end
+
+    if m < 2*k + 1
+        throw(DomainError(m, "m must be >= 2*k + 1"))
+    end
+
+    D = zeros(m+2,m+1)
+    if k == 2 
+        for i = 2:(m+1)
+            D[i,(i-1):i] = [-1 1]
+        end
+    elseif k == 4
+        A = [-11/12 17/24 3/8 -5/24 1/24]
+        D[2,1:5] = A
+        D[m+1, (m-3):(m+1)] = -reverse(A)
+        for i = 3:m
+            D[i,(i-2):(i+1)] = [1/24 -9/8 9/8 -1/24]
+        end
+    elseif k == 6
+        A = [-1627/1920  211/640  59/48  -235/192 91/128 -443/1920 31/960;
+                31/960  -687/640 129/128   19/192 -3/32    21/640  -3/640]
+        D[2:3,1:7] = A
+        D[m:(m+1), (m-5):(m+1)] = -rot180(A)
+        for i = 4:(m-1)
+            D[i,(i-3):(i+2)] = [-3/640 25/384 -75/64 75/64 -25/384 3/640]
+        end
+    elseif k == 8
+        A = [-1423/1792     -491/7168   7753/3072 -18509/5120  3535/1024 -2279/1024  953/1024 -1637/7168  2689/107520;
+              2689/107520 -36527/35840  4259/5120   6497/15360 -475/1024  1541/5120 -639/5120  1087/35840  -59/17920;
+               -59/17920    1175/21504 -1165/1024   1135/1024    25/3072  -251/5120   25/1024   -45/7168     5/7168]
+            D[2:4, 1:9] = A;
+        D[2:4,1:9] = A
+        D[(m-1):(m+1), (m-7):(m+1)] = -rot180(A)
+        for i = 5:(m-2)
+            D[i,(i-4):(i+3)] = [5/7168 -49/5120 245/3072 -1225/1024 1225/1024 -245/3072 49/5120 -5/7168]
+        end
+    end
+    D = (1/dx)*D;
+end
\ No newline at end of file
diff --git a/julia/MOLE.jl/src/gradient.jl b/julia/MOLE.jl/src/gradient.jl
new file mode 100644
index 00000000..51b66067
--- /dev/null
+++ b/julia/MOLE.jl/src/gradient.jl
@@ -0,0 +1,61 @@
+"""
+    grad(k, m, dx)
+
+Returns a m+1 by m+2 one-dimensional mimetic gradient operator.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells
+- `dx`: Step size
+"""
+function grad(k::Int,m::Int,dx)
+    if k < 2
+        throw(DomainError(k, "k must be >= 2"))
+    end
+
+    if k % 2 != 0
+        throw(DomainError(k, "k must be an positive even integer"))
+    end
+
+    if m < 2*k
+        throw(DomainError(m, "m must be >= 2*k"))
+    end
+
+    G = zeros(m+1,m+2)
+    if k == 2
+        A = [-8/3 3 -1/3]
+        G[1,1:3] = A #[-8/3 3 -1/3]
+        G[m+1,m:(m+2)] = -reverse(A) #[1/3, -3, 8/3]
+        for i = 2:m
+            G[i,i:(i+1)] = [-1 1]
+        end
+    elseif k == 4
+        A = [-352/105  35/8  -35/24  21/40  -5/56; 
+              16/105  -31/24  29/24  -3/40   1/168]
+        G[1:2,1:5] = A
+        G[m:(m+1), (m-2):(m+2)] = -rot180(A)
+        for i = 3:(m-1)
+            G[i,(i-1):(i+2)] = [1/24 -9/8 9/8 -1/24]
+        end
+    elseif k == 6
+        A = [-13016/3465  693/128  -385/128  693/320  -495/448  385/1152  -63/1408;
+                496/3465 -811/640   449/384  -29/960   -11/448   13/1152  -37/21120;
+                 -8/385   179/1920 -153/128  381/320  -101/1344   1/128    -3/7040];
+        G[1:3,1:7] = A
+        G[(m-1):(m+1), (m-4):(m+2)] = -rot180(A)
+        for i = 4:(m-2)
+            G[i,(i-2):(i+3)] = [-3/640 25/384 -75/64 75/64 -25/384 3/640]
+        end
+    elseif k == 8
+        A = [-182144/45045     6435/1024    -5005/1024  27027/5120  -32175/7168  25025/9216  -12285/11264  3465/13312   -143/5120;
+               86048/675675 -131093/107520  49087/46080 10973/76800  -4597/21504  4019/27648 -10331/168960 2983/199680 -2621/1612800;
+               -3776/225225    8707/107520 -17947/15360 29319/25600   -533/21504  -263/9216     903/56320  -283/66560    257/537600;
+                  32/9009      -543/35840     265/3072  -1233/1024    8625/7168   -775/9216     639/56320  -15/13312       1/21504]
+        G[1:4,1:9] = A
+        G[(m-2):(m+1), (m-6):(m+2)] = -rot180(A)
+        for i = 5:(m-3)
+            G[i,(i-3):(i+4)] = [5/7168 -49/5120 245/3072 -1225/1024 1225/1024 -245/3072 49/5120 -5/7168]
+        end
+    end
+    G = (1/dx)*G;
+end
\ No newline at end of file
diff --git a/julia/MOLE.jl/src/laplacian.jl b/julia/MOLE.jl/src/laplacian.jl
new file mode 100644
index 00000000..ea625e22
--- /dev/null
+++ b/julia/MOLE.jl/src/laplacian.jl
@@ -0,0 +1,16 @@
+"""
+    lap(k, m, dx)
+
+Returns a m+2 by m+2 one-dimensional mimetic laplacian operator.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells
+- `dx`: Step size
+"""
+function lap(k::Int, m::Int, dx)
+    D = div(k,m,dx)
+    G = grad(k,m,dx)
+
+    L = D*G;
+end
\ No newline at end of file
diff --git a/julia/MOLE.jl/src/robinBC.jl b/julia/MOLE.jl/src/robinBC.jl
new file mode 100644
index 00000000..a0283278
--- /dev/null
+++ b/julia/MOLE.jl/src/robinBC.jl
@@ -0,0 +1,25 @@
+"""
+    robinBC(k, m, dx, a, b)
+
+Returns a m+2 by m+2 one-dimensional mimetic boundary operator that imposes a boundary condition of Robin's type.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells
+- `dx`: Step size
+- `a`: Dirichlet Coefficient
+- `b`: Neumann Coefficient
+"""
+function robinBC(k::Int, m::Int, dx, a, b)
+    A = zeros(m+2,m+2)
+    A[1,1] = a
+    A[m+2,m+2] = a
+
+    B = zeros(m+2,m+1)
+    B[1,1] = -b
+    B[m+2,m+1] = b
+
+    G = grad(k,m,dx)
+
+    BC = A + B*G;
+end
\ No newline at end of file
diff --git a/julia/MOLE.jl/test/Project.toml b/julia/MOLE.jl/test/Project.toml
new file mode 100644
index 00000000..a2cd2f8d
--- /dev/null
+++ b/julia/MOLE.jl/test/Project.toml
@@ -0,0 +1,3 @@
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/julia/MOLE.jl/test/runtests.jl b/julia/MOLE.jl/test/runtests.jl
new file mode 100644
index 00000000..86b8000a
--- /dev/null
+++ b/julia/MOLE.jl/test/runtests.jl
@@ -0,0 +1,26 @@
+using MOLE
+using Test, LinearAlgebra
+
+#tol = 1e-10
+
+@testset "Testing MOLE operators" begin
+    @testset "Testing 1D divergence" begin
+        include("testDivergence.jl")
+    end
+
+    @testset "Testing 1D gradient" begin
+        include("testGradient.jl")
+    end
+
+    @testset "Testing 1D laplacian" begin
+        include("testLaplacian.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
diff --git a/julia/MOLE.jl/test/testDivergence.jl b/julia/MOLE.jl/test/testDivergence.jl
new file mode 100644
index 00000000..9a0dd6d6
--- /dev/null
+++ b/julia/MOLE.jl/test/testDivergence.jl
@@ -0,0 +1,9 @@
+tol = 1e-10
+
+@testset "Testing divergence for order k=$k" for k=2:2:8
+    m = 2*k+1
+    D = MOLE.div(k, m, 1/m)
+    field = ones(m+1,1)
+    sol = D*field
+    @test norm(sol) < tol
+end
\ No newline at end of file
diff --git a/julia/MOLE.jl/test/testGradient.jl b/julia/MOLE.jl/test/testGradient.jl
new file mode 100644
index 00000000..3cc0ffd4
--- /dev/null
+++ b/julia/MOLE.jl/test/testGradient.jl
@@ -0,0 +1,9 @@
+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)
+    field = ones(m+2, 1)
+    sol = G*field
+    @test norm(sol) < tol
+end;
\ No newline at end of file
diff --git a/julia/MOLE.jl/test/testLaplacian.jl b/julia/MOLE.jl/test/testLaplacian.jl
new file mode 100644
index 00000000..8103589c
--- /dev/null
+++ b/julia/MOLE.jl/test/testLaplacian.jl
@@ -0,0 +1,9 @@
+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)
+    field = ones(m+2, 1)
+    sol = L*field
+    @test norm(sol) < tol
+end;
\ No newline at end of file