diff --git a/julia/MOLE.jl/docs/make.jl b/julia/MOLE.jl/docs/make.jl
index 5c8325a2..37473b12 100644
--- a/julia/MOLE.jl/docs/make.jl
+++ b/julia/MOLE.jl/docs/make.jl
@@ -1,5 +1,5 @@
 push!(LOAD_PATH,"../src/")
-using Documenter, MOLE
+using Documenter, MOLE, SparseArrays
 
 makedocs(
     sitename="MOLE.jl Docs",
diff --git a/julia/MOLE.jl/docs/src/index.md b/julia/MOLE.jl/docs/src/index.md
index afab4749..7f87d49b 100644
--- a/julia/MOLE.jl/docs/src/index.md
+++ b/julia/MOLE.jl/docs/src/index.md
@@ -45,6 +45,7 @@ julia --project=. path/to/myScript.jl
 To run the unit tests, first enter the Julia REPL as in the above section (that is, by running the command `julia --project=.` from the directory `mole/julia/MOLE.jl`). Next, enter the `pkg` mode by pressing `]`, then type the command `test`. The results of the unit tests should be displayed to your console.
 
 ## Building the documentation
+
 MOLE.jl uses [Documenter.jl](https://documenter.juliadocs.org/stable/) to build its Julia implementation documentation. From the `mole/julia/MOLE.jl` directory, navigate to the `docs/` directory, with 
 
 ```sh
@@ -62,6 +63,7 @@ include("make.jl")
 ```
 
 ### From the command line
+
 To build the documentation from the `docs/` directory, use the following commands to instatiate and precompile the documentation environment:
 
 ```sh
@@ -79,6 +81,7 @@ julia --project=. make.jl
 ```
 
 ### Preview the documentation
+
 Once you have built the documentation (either from the REPL or the command line), you can inspect the documentation in `docs/build` with the `index.html` file as the homepage.
 
 ## Functions
@@ -86,15 +89,27 @@ Once you have built the documentation (either from the REPL or the command line)
 ### Operators
 
 ```@docs
-MOLE.div(k::Int,m::Int,dx)
-MOLE.grad(k::Int,m::Int,dx)
-MOLE.lap(k::Int,m::Int,dx)
+MOLE.Operators.div(k::Int, m::Int, dx; dc::NTuple{2,T}, nc::NTuple{2,T})
+MOLE.Operators.div(k::Int, ticks::AbstractVector)
+MOLE.Operators.div(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4,T}, nc::NTuple{4,T})
+MOLE.Operators.div(k::Int, xticks::AbstractVector, yticks::AbstractVector)
+MOLE.Operators.grad(k::Int, m::Int, dx; dc::NTuple{2,T}, nc::NTuple{2,T})
+MOLE.Operators.grad(k::Int, ticks::AbstractVector)
+MOLE.Operators.grad(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4,T}, nc::NTuple{4,T})
+MOLE.Operators.grad(k::Int, xticks::AbstractVector, yticks::AbstractVector)
+MOLE.Operators.lap(k::Int, m::Int, dx; dc::NTuple{2,T}, nc::NTuple{2,T})
+MOLE.Operators.lap(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4,T}, nc::NTuple{4,T})
 ```
 
 ### Utilities
 
 ```@docs
-MOLE.robinBC(k::Int, m::Int, dx, a, b)
+MOLE.BCs.robinBC(k::Int, m::Int, dx, a, b)
+MOLE.BCs.robinBC(k::Int, m::Int, dx, n::Int, dy, a, b)
+MOLE.BCs.ScalarBC1D{T}
+MOLE.BCs.ScalarBC2D{T}
+MOLE.BCs.addScalarBC!(A::SparseMatrixCSC, b::AbstractVector, bc::MOLE.BCs.ScalarBC1D, k::Integer, m::Integer, dx)
+MOLE.BCs.addScalarBC!(A::SparseMatrixCSC, b::AbstractVector, bc::MOLE.BCs.ScalarBC2D, k::Integer, m::Integer, dx, n::Integer, dy)
 ```
 
 ## Examples
@@ -104,5 +119,11 @@ The MOLE library contains examples demonstrating how to use the operators, in a
 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
+  - 1D Examples
+    - ```elliptic1D```: A script that solves the 1D Poisson's equation with Robin boundary conditions using mimetic operators.
+  - 2D Examples
+    - ```elliptic2DXDirichletYDirichlet```: A script that solves the 2D Laplace equation, $\nabla^2 u = 0$, with Dirichlet boundary conditions in $x$ and $y$ using mimetic operators.
+    - ```elliptic2DXPerYDirichlet```: A script that solves the 2D Laplace equation, $\nabla^2 u = 0$, with periodic bonudary conditions in $x$ and Dirichlet boundary conditions in $y$ using mimetic operators.
+- Parabolic Problems
+  - 2D Examples
+    - ```parabolic2D```: A script that solves the 2D heat equation, $u_t = \nu \nabla^2 u$, with Dirichlet boundary conditions in $x$ and $y$ using mimetic operators.
diff --git a/julia/MOLE.jl/examples/elliptic2DXDirichletYDirichlet.jl b/julia/MOLE.jl/examples/elliptic2DXDirichletYDirichlet.jl
new file mode 100644
index 00000000..04cb688a
--- /dev/null
+++ b/julia/MOLE.jl/examples/elliptic2DXDirichletYDirichlet.jl
@@ -0,0 +1,56 @@
+using LinearAlgebra
+using SparseArrays
+using Plots
+
+import MOLE: Operators, BCs
+
+# Parameters
+k = 2       # Order of accuracy
+m = 99      # Number of cells in x-direction
+n = m + 2   # Number of cells in y-direction
+dx = pi / m # Step size in x-direction
+dy = pi / n # Step size in y-direction
+
+# Grid
+xc = [0; (dx / 2.0) : dx : (pi - dx / 2.0); pi]
+yc = [0; (dy / 2.0) : dy : (pi - dy / 2.0); pi]
+X = ones(1, n + 2) .* xc
+Y = yc' .* ones(m + 2, 1)
+
+# Exact Solution
+ue = exp.(X) .* cos.(Y)
+
+# Boundary Conditions
+dc = (1.0, 1.0, 1.0, 1.0)
+nc = (0.0, 0.0, 0.0, 0.0)
+v = (vec(ue[1,2:end-1]'), vec(ue[end,2:end-1]'), vec(ue[:,1]), vec(ue[:, end]))
+bc = BCs.ScalarBC2D(dc, nc, v)
+
+# Operator
+A = - Operators.lap(k, m, dx, n, dy, dc=dc, nc=nc)
+
+# RHS
+b = zeros(m + 2, n + 2)
+b = vec(b)
+
+# Apply BCs
+A0, b0 = BCs.addScalarBC!(sparse(A), b, bc, k, m, dx, n, dy)
+
+# Approximate Solution
+ua = A0 \ b0
+ua = Matrix(reshape(ua, m + 2, n + 2))
+
+pa = heatmap(xc, yc, ua, title = "Approximate Solution", xlabel = "X", ylabel = "Y", colorbar_title = "u(x,y)", aspect_ratio = :equal)
+display(pa)
+println("Press Enter to close the plot and open the next.")
+readline()
+
+pe = heatmap(xc, yc, ue, title = "Exact Solution", xlabel = "X", ylabel = "Y", colorbar_title = "u(x,y)", aspect_ratio = :equal)
+display(pe)
+println("Press Enter to close the plot and exit.")
+readline()
+
+max_err = maximum(abs, ue - ua)
+println("Maximum error: $max_err")
+rel_err = 100 * max_err ./ (maximum(ue) - minimum(ue))
+println("Relative error: $rel_err")
\ No newline at end of file
diff --git a/julia/MOLE.jl/examples/elliptic2DXPerYDirichlet.jl b/julia/MOLE.jl/examples/elliptic2DXPerYDirichlet.jl
new file mode 100644
index 00000000..26e45547
--- /dev/null
+++ b/julia/MOLE.jl/examples/elliptic2DXPerYDirichlet.jl
@@ -0,0 +1,51 @@
+using LinearAlgebra
+using SparseArrays
+using Plots
+
+import MOLE: Operators, BCs
+
+# Parameters
+k = 2        # Order of accuracy
+m = 20       # Number of cells in x-direction
+n = m + 1    # Number of cells in y-direction
+dx = 1.0 / m # Step size in x-direction
+dy = 1.0 / n # Step size in y-direction
+
+# Grid
+xc = (dx / 2.0) : dx : (1.0 - dx / 2.0)
+yc = [0; (dy / 2.0) : dy : (1.0 - dy / 2.0); 1.0]
+X = ones(1, n + 2) .* xc
+Y = yc' .* ones(m, 1)
+
+# Exact Solution
+ue = Y .* (1.0 .- Y) .* sin.(2.0 * pi .* X)
+
+# Boundary Conditions
+dc = (0.0, 0.0, 1.0, 1.0)
+nc = (0.0, 0.0, 0.0, 0.0)
+v = ([0.0], [0.0], vec(zeros(m, 1)), vec(zeros(m, 1)))
+bc = BCs.ScalarBC2D(dc, nc, v)
+
+# Operator
+A = - Operators.lap(k, m, dx, n, dy, dc=dc, nc=nc)
+
+# RHS
+b = 2.0 * sin.(2.0 * pi .* X) .* (1.0 .+ 2.0 * pi^2 .* Y .* (1.0 .- Y))
+b = vec(b)
+
+# Apply BCs
+A0, b0 = BCs.addScalarBC!(sparse(A), b, bc, k, m, dx, n, dy)
+
+# Approximate Solution
+ua = A0 \ b0
+ua = Matrix(reshape(ua, m, n + 2))
+
+pa = heatmap(xc, yc, ua, title = "Approximate Solution", xlabel = "X", ylabel = "Y", colorbar_title = "u(x,y)", aspect_ratio = :equal)
+display(pa)
+println("Press Enter to close the plot and open the next.")
+readline()
+
+pe = heatmap(xc, yc, ue, title = "Exact Solution", xlabel = "X", ylabel = "Y", colorbar_title = "u(x,y)", aspect_ratio = :equal)
+display(pe)
+println("Press Enter to close the plot and exit.")
+readline()
\ No newline at end of file
diff --git a/julia/MOLE.jl/examples/parabolic2D.jl b/julia/MOLE.jl/examples/parabolic2D.jl
new file mode 100644
index 00000000..9c2e6fe9
--- /dev/null
+++ b/julia/MOLE.jl/examples/parabolic2D.jl
@@ -0,0 +1,74 @@
+using LinearAlgebra
+using SparseArrays
+using Plots
+
+import MOLE: Operators, BCs
+
+# Parameters
+k = 2               # Order of accuracy
+nu = 0.1            # Diffusion coefficient
+m = 40              # Number of cells in x-direction
+n = 50              # Number of cells in y-direction
+dx = 2.0 / m        # Step size in x-direction
+dy = 2.0 / n        # Step size in y-direction
+t = 3.0             # Final time
+method = "implicit" # Method of Euler time integration
+dt = 0.0            # Step size in time
+if method == "explicit"; dt = 0.001; else; dt = 0.01; end
+
+# Grid
+xc = [0; (dx / 2.0) : dx : (2.0 - dx / 2.0); 2.0]
+yc = [0; (dy / 2.0) : dy : (2.0 - dy / 2.0); 2.0]
+
+# Initial Conditions
+u = zeros(m + 2, n + 2)
+for i = 1:m
+    for j = 1:n
+        if ((1.0 ≤ yc[j]) && (yc[j] ≤ 1.5) && (1.0 ≤ xc[i]) && (xc[i] ≤ 1.5))
+            u[i, j] = 2.0
+        end
+    end
+end
+u = Matrix(reshape(u, :, 1))
+
+# Boundary Conditions
+dc = (1.0, 1.0, 1.0, 1.0)
+nc = (0.0, 0.0, 0.0, 0.0)
+v = (vec(zeros(n, 1)), vec(zeros(n, 1)), vec(zeros(m + 2, 1)), vec(zeros(m + 2, 1)))
+bc = BCs.ScalarBC2D(dc, nc, v)
+
+# Operator
+L = Operators.lap(k, m, dx, n, dy, dc=dc, nc=nc)
+if method == "explicit"
+    L = nu * dt .* sparse(L) .+ sparse(Matrix(I, size(L)))
+else
+    L = sparse(Matrix(I, size(L))) .- nu * dt .* sparse(L)
+end
+
+for it = 0:(t/dt)
+    heatmap(xc, yc, u,
+        size = (800, 600),
+        xlabel = "x",
+        ylabel = "y",
+        title = "$method\n2-D Diffusion with ν = $nu -- time (t) = $(it*dt)",
+        xlims = (0, 2),
+        ylims = (0, 2),
+        clims = (0, 2),
+        colorbar = true,
+        colorbar_title = "u(x,y)"
+    )
+    display(current())
+
+    global u = reshape(u', :, 1)
+
+    if method == "explicit"
+        global u = L * u
+    else
+        global u = L \ u
+    end
+
+    global u = reshape(u, n + 2, m + 2)'
+end
+
+println("Press Enter to close the plot and exit.")
+readline()
\ No newline at end of file
diff --git a/julia/MOLE.jl/src/BCs/robinBC.jl b/julia/MOLE.jl/src/BCs/robinBC.jl
index 0fbd5e57..9bd266cd 100644
--- a/julia/MOLE.jl/src/BCs/robinBC.jl
+++ b/julia/MOLE.jl/src/BCs/robinBC.jl
@@ -30,4 +30,45 @@ function robinBC(k::Int, m::Int, dx, a, b)
     G = Operators.grad(k,m,dx)
 
     BC = A + B*G;
+end
+
+"""
+    robinBC2D(k, m, dx, n, dy, a, b)
+
+Returns a two-dimensional mimetic boundary condition operator that imposes a boundary condition of Robin's type
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells in x-direction
+- `dx`: Step size in x-direction
+- `n::Int`: Number of cells in y-direction
+- `dy`: Step size in y-direction
+- `a`: Dirichlet Coefficient
+- `b`: Neumann Coefficient
+"""
+function robinBC2D(k::Int, m::Int, dx, n::Int, dy, a, b)
+
+    Bm = robinBC(k, m, dx, a, b)
+    Bn = robinBC(k, n, dy, a, b)
+
+    Im = Matrix(I, m + 2, m + 2)
+
+    In = Matrix(I, n + 2, n + 2)
+    In[1, 1] = 0
+    In[end, end] = 0
+
+    BC1 = kron(In, Bm)
+    BC2 = kron(Bn, Im)
+
+    BC = BC1 + BC2
+
+end
+
+"""
+    robinBC(k, m, dx, n, dy, a, b)
+
+Alias of robinBC2D
+"""
+function robinBC(k::Int, m::Int, dx, n::Int, dy, a, b)
+    return robinBC2D(k, m, dx, n, dy, a, b);
 end
\ No newline at end of file
diff --git a/julia/MOLE.jl/src/BCs/scalarBC.jl b/julia/MOLE.jl/src/BCs/scalarBC.jl
index c28ecf80..1bf6ce60 100644
--- a/julia/MOLE.jl/src/BCs/scalarBC.jl
+++ b/julia/MOLE.jl/src/BCs/scalarBC.jl
@@ -4,6 +4,7 @@
     See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details.
 =#
 
+using LinearAlgebra
 import ..Operators: Operators, grad
 
 """
@@ -28,6 +29,22 @@ struct ScalarBC1D{T} <: AbstractScalarBC{1}
     v::NTuple{2,T}
 end
 
+
+"""
+    Concrete scalar BC description for 2D.
+
+    Fields mirror the MATLAB function:
+    - dc: Dirichlet coefficients (left, right, bottom, top)
+    - nc: Neumann/Robin coefficients (left, right, bottom, top)
+    - v: prescribed boundary value g (left, right, bottom, top)
+"""
+struct ScalarBC2D{T} <: AbstractScalarBC{2}
+    dc::NTuple{4,T}
+    nc::NTuple{4,T}
+    v::NTuple{4,AbstractVector{T}}
+end
+
+
 # -------------------------
 # 1D implementation
 # -------------------------
@@ -104,6 +121,154 @@ function addScalarBC!(A::SparseMatrixCSC, b::AbstractVector, bc::ScalarBC1D{T},
 end
 
 
+# -------------------------
+# 2D implementation
+# -------------------------
+
+"""
+    Internal helper: build LHS contributions for 2D BC.
+    Returns (Al, Ar, Ab, At).
+"""
+function _scalarbc2D_lhs(k::Integer, m::Integer, dx, n::Integer, dy, dc::NTuple{4,T}, nc::NTuple{4,T}) where {T}
+
+    Abcl, Abcr, Abcb, Abct = 0, 0, 0, 0 # periodic case
+
+    hasbclr = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T))
+    hasbcbt = (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T))
+
+    if hasbclr
+
+        Abcl0, Abcr0 = _scalarbc1d_lhs(k, m, dx, dc[1:2], nc[1:2])
+
+        if !hasbcbt
+            In = sparse(Matrix(I, n, n))
+        else
+            In = Matrix(I, n + 2, n + 2)
+            In[1, 1] = 0
+            In[end, end] = 0
+            In = sparse(In)
+        end
+
+        # Left and right edges
+        Abcl = kron(In, Abcl0)
+        Abcr = kron(In, Abcr0)
+
+    end
+
+    if hasbcbt
+
+        Abcb0, Abct0 = _scalarbc1d_lhs(k, n, dy, dc[3:4], nc[3:4])
+
+        if !hasbclr
+            Im = sparse(Matrix(I, m, m))
+        else
+            Im = sparse(Matrix(I, m + 2, m + 2))
+        end
+
+        Abcb = kron(Abcb0, Im)
+        Abct = kron(Abct0, Im)
+
+    end
+
+    return Abcl, Abcr, Abcb, Abct
+
+end
+
+"""
+    Internal helper: overwrite RHS at boundary indices
+"""
+@inline function _scalarbc2d_rhs(b, dc::NTuple{4,T}, nc::NTuple{4,T}, v::NTuple{4,AbstractVector{T}}, rl::AbstractVector{Int}, rr::AbstractVector{Int}, rb::AbstractVector{Int}, rt::AbstractVector{Int}) where {T}
+
+    hasbclr = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T))
+    hasbcbt = (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T))
+
+    if hasbclr
+        b[rl] = v[1]
+        b[rr] = v[2]
+    end
+
+    if hasbcbt
+        b[rb] = v[3]
+        b[rt] = v[4]
+    end
+
+    return b
+end
+
+"""
+    2D BC applicator. Mirrors MATLAB addScalarBC2D.
+    Signature keeps the discretization params (`k,m,dx,n,dy`) separate from `bc`.
+"""
+function addScalarBC!(A::SparseMatrixCSC, b::AbstractVector, bc::ScalarBC2D{T},
+                      k::Integer, m::Integer, dx, n::Integer, dy) where {T}
+
+    dc, nc, v = bc.dc, bc.nc, bc.v
+
+    hasbclr = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T))
+    hasbcbt = (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T))
+
+    rl, rr, rb, rt = [0], [0], [0], [0] # periodic case
+
+    Abcl, Abcr, Abcb, Abct = _scalarbc2D_lhs(k, m, dx, n, dy, dc, nc)
+
+    if hasbclr
+
+        rl, _, _ = findnz(Abcl)
+        rr, _, _ = findnz(Abcr)
+        rl = unique(rl)
+        rr = unique(rr)
+
+        # remove rows of A associated to boundary
+        Abc1 = Abcl .+ Abcr
+        rowsbc1, _, _ = findnz(Abc1)
+        # For some reason, the below code doesn't work:
+            # rows1, cols1, s1 = findnz(A[rowsbc1, :])
+            # A = A .- sparse(rows1, cols1, s1, size(A, 1), size(A, 2))
+        A[rowsbc1, :] = spzeros(size(A[rowsbc1, :])) # This does though
+        # Update matrix A with boundary information
+        A = A .+ Abc1
+        # Remove b entries associated to bcs
+        b[rowsbc1] .= 0
+
+    end
+
+    if hasbcbt
+        
+        rb, _, _ = findnz(Abcb)
+        rt, _, _ = findnz(Abct)
+        rb = unique(rb)
+        rt = unique(rt)
+
+        # Remove rows of A associated to boundary
+        Abc2 = Abct .+ Abcb
+        rowsbc2, _, _ = findnz(Abc2)
+        # rows2, cols2, s2 = findnz(A[rowsbc2, :])
+        # A = A .- sparse(rows2, cols2, s2, size(A, 1), size(A, 2))
+        A[rowsbc2, :] = spzeros(size(A[rowsbc2, :]))
+        # Update matriz A with boundary information
+        A = A .+ Abc2
+        # Remove b entries associated to bcs
+        b[rowsbc2] .= 0
+
+    end
+
+    if hasbclr || hasbcbt
+        b = _scalarbc2d_rhs(b, dc, nc, v, rl, rr, rb, rt)
+    end
+
+    return A, b
+end
+
+
+"""
+    Alias for `addScalarBC!(A, b, bc, k, m, dx, n, dy)`
+"""
+function addScalarBC2D!(A::SparseMatrixCSC, b::AbstractVector{T}, bc::ScalarBC2D{T},
+                      k::Integer, m::Integer, dx, n::Integer, dy) where {T}
+    return addScalarBC!(A, b, bc, k, m, dx, n, dy)
+end
+
+
 # -------------------------
 # Construction helpers
 # -------------------------
diff --git a/julia/MOLE.jl/src/Operators/divergence.jl b/julia/MOLE.jl/src/Operators/divergence.jl
index 238cd347..7d5af46f 100644
--- a/julia/MOLE.jl/src/Operators/divergence.jl
+++ b/julia/MOLE.jl/src/Operators/divergence.jl
@@ -4,8 +4,50 @@
     See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details.
 =#
 
+# ------------------------
+# 1-D Divergence Operators
+# ------------------------
+
 """
-    div(k, m, dx)
+    div(k, m, dx; dc, nc)
+
+Returns a one-dimensional mimetic divergence operator. Default is non periodic.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells
+- `dx::T`: Step size
+- `dc::NTuple{2,T}`: Dirichlet coefficients of left and right boundaries (optional)
+- `nc::NTuple{2,T}`: Neumann coefficients of left and right boundaries (optional)
+"""
+function div(k::Int, m::Int, dx; dc::NTuple{2,T} = (1.0, 1.0), nc::NTuple{2,T} = (1.0, 1.0)) where {T}
+
+    hasbc = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T))
+    if hasbc
+        D = divNonPeriodic(k, m, dx)
+        return D
+    else
+        D = divPeriodic(k, m, dx)
+        return D
+    end
+
+end
+
+"""
+    div(k, ticks)
+
+Returns a m + 2 by m + 1 non-uniform mimetic divergence operator
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `ticks::AbstractArray`: Edges' ticks e.g. [0 0.1 0.15 0.2 0.3 0.4 0.45]
+"""
+function div(k::Int, ticks::AbstractArray)
+    return divNonUniform(k, ticks)
+end
+
+"""
+    divNonPeriodic(k, m, dx)
 
 Returns a m+2 by m+1 one-dimensional mimetic divergence operator.
 
@@ -14,9 +56,9 @@ Returns a m+2 by m+1 one-dimensional mimetic divergence operator.
 - `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"))
+function divNonPeriodic(k::Int, m::Int, dx)
+    if k < 2 || k > 8
+        throw(DomainError(k, "k must be >= 2 and <= 8"))
     end
 
     if k % 2 != 0
@@ -59,4 +101,144 @@ function div(k::Int,m::Int,dx)
         end
     end
     D = (1/dx)*D;
+end
+
+
+"""
+    divPeriodic(k, m, dx)
+
+Returns a m by m periodic mimetic divergence operator.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells
+- `dx`: Step size
+"""
+function divPeriodic(k::Int, m::Int, dx)
+
+    D = - gradPeriodic(k, m, dx)';
+    
+end
+
+
+"""
+    divNonUniform(k, ticks)
+
+Returns a m + 2 by m + 1 non-uniform mimetic divergence operator
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `ticks::AbstractVector`: Edges' ticks e.g. [0 0.1 0.15 0.2 0.3 0.4 0.45]
+"""
+function divNonUniform(k::Int, ticks::AbstractVector)
+
+    D = div(k, length(ticks) - 1, 1)
+
+    m, _ = size(D)
+
+    if size(ticks, 1) == 1
+        J = diagm((D * ticks').^-1)
+    else
+        J = diagm((D * ticks).^-1)
+    end
+
+    D = J * D
+    D[1, :] .= 0
+    D[end, :] .= 0
+    return D;
+
+end
+
+# ------------------------
+# 2-D Divergence Operators
+# ------------------------
+
+"""
+    div(k, m, dx, n, dy; dc, nc)
+
+Returns a two-dimensional mimetic divergence operator. Default is non periodic.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells in x-direction
+- `dx::T`: Step size in x-direction
+- `n::Int`: Number of cells in y-direction
+- `dy::T`: Step size in y-direction
+- `dc::NTuple{4,T}`: Dirichlet coefficients of left, right, bottom, and top boundaries (optional)
+- `nc::NTuple{4,T}`: Neumann coefficients of left, right, bottom, and top boundaries (optional)
+"""
+function div(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4,T} = (1.0, 1.0, 1.0, 1.0), nc::NTuple{4,T} = (1.0, 1.0, 1.0, 1.0)) where {T}
+
+    hasbclr = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T))
+    hasbcbt = (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T))
+
+    if hasbclr
+        Dx = divNonPeriodic(k, m, dx)
+        Im = Matrix(I, m + 2, m + 2)
+        Im = Im[:, 2:end-1]
+    else
+        Dx = divPeriodic(k, m, dx)
+        Im = Matrix(I, m, m)
+    end
+
+    if hasbcbt
+        Dy = divNonPeriodic(k, n, dy)
+        In = Matrix(I, n + 2, n + 2)
+        In = In[:, 2:end-1]
+    else
+        Dy = divPeriodic(k, n, dy)
+        In = Matrix(I, n, n)
+    end
+
+    Sx = kron(In, Dx)
+    Sy = kron(Dy, Im)
+
+    D = [Sx Sy];
+
+end
+
+"""
+    div(k, xticks, yticks)
+
+Returns a two-dimensional non-uniform mimetic divergence operator.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `xticks::AbstractVector`: Edges' ticks (x-axis)
+- `yticks::AbstractVector`: Edges' ticks (y-axis)
+"""
+function div(k::Int, xticks::AbstractVector, yticks::AbstractVector)
+    return div2DNonUniform(k, xticks, yticks)
+end
+
+
+"""
+    div2DNonUniform(k, xticks, yticks)
+
+Returns a two-dimensional non-uniform mimetic divergence operator.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `xticks::AbstractVector`: Edges' ticks (x-axis)
+- `yticks::AbstractVector`: Edges' ticks (y-axis)
+"""
+function div2DNonUniform(k::Int, xticks::AbstractVector, yticks::AbstractVector)
+
+    Dx = divNonUniform(k, xticks)
+    Dy = divNonUniform(k, yticks)
+
+    m = size(Dx, 1) # Really m + 2, but makes for simpler augmented identity matrix constuction
+    n = size(Dy, 1) # Really n + 2, but makes for simpler augmented identity matrix constuction
+
+    Im = Matrix(I, m, m)
+    In = Matrix(I, n, n)
+
+    Im = Im[:, 2:end-1]
+    In = In[:, 2:end-1]
+
+    Sx = kron(In, Dx)
+    Sy = kron(Dy, Im)
+
+    D = [Sx Sy]
+
 end
\ No newline at end of file
diff --git a/julia/MOLE.jl/src/Operators/gradient.jl b/julia/MOLE.jl/src/Operators/gradient.jl
index 24b0af3f..0f6a2a21 100644
--- a/julia/MOLE.jl/src/Operators/gradient.jl
+++ b/julia/MOLE.jl/src/Operators/gradient.jl
@@ -4,6 +4,52 @@
     See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details.
 =#
 
+using LinearAlgebra
+
+#-----------------------
+# 1-D Gradient Operators
+#-----------------------
+
+"""
+    grad(k, m, dx; dc, nc)
+
+Returns a one-dimensional mimetic gradient operator. Default is non periodic.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells
+- `dx::T`: Step size
+- `dc::NTuple{2,T}`: Dirichlet coefficients of the left and right boundaries (optional)
+- `nc::NTuple{2,T}`: Neumann coefficients of the left and right boundaries (optional)
+"""
+function grad(k::Int, m::Int, dx; dc::NTuple{2,T} = (1.0, 1.0), nc::NTuple{2,T} = (1.0, 1.0)) where {T}
+
+    hasbc = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T))
+    if hasbc
+        G = gradNonPeriodic(k, m, dx)
+        return G
+    else
+        G = gradPeriodic(k, m, dx)
+        return G
+    end
+
+end
+
+
+"""
+    grad(k, ticks)
+
+Returns a m + 1 by m + 2 one-dimensional non-uniform mimetic gradient operator
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `ticks::AbstractVector`: Centers' ticks e.g. (0 0.5 1 3 5 7 9 9.5 10) (including the boundaries!)
+"""
+function grad(k::Int, ticks::AbstractVector)
+    return gradNonUniform(k, ticks)
+end
+
+
 """
     grad(k, m, dx)
 
@@ -14,9 +60,9 @@ Returns a m+1 by m+2 one-dimensional mimetic gradient operator.
 - `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"))
+function gradNonPeriodic(k::Int, m::Int, dx)
+    if k < 2 || k > 8
+        throw(DomainError(k, "k must be >= 2 and <= 8"))
     end
 
     if k % 2 != 0
@@ -64,4 +110,238 @@ function grad(k::Int,m::Int,dx)
         end
     end
     G = (1/dx)*G;
+end
+
+
+"""
+    gradPeriodic(k, m, dx)
+
+Returns a m by m periodic mimetic gradient operator
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells in x-direction
+- `dx`: Step size in x-direction
+"""
+function gradPeriodic(k::Int, m::Int, dx)
+    
+   if k < 2 || k > 8
+        throw(DomainError(k, "k must be >= 2 and <= 8"))
+    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
+    
+    V = zeros(1, m)
+    idx = fill(-1, m, m)
+    idx[:, 1] = 1:m
+    idx = cumsum(idx, dims=2)
+    idx = mod.(idx .+ m, m) .+ 1
+
+    if k == 2
+
+        V[1, 2:3] = [1, -1]
+
+    elseif k == 4
+
+        V[1, 1:4] = [-1/24, 9/8, -9/8, 1/24]
+
+    elseif k == 6
+
+        V[1, 1:5] = [-25/384, 75/64, -75/64, 25/384, -3/640]
+        V[1, end] = 3/640
+
+    elseif k == 8
+
+        V[1, 1:6] = [-245/3072, 1225/1024, -1225/1024, 245/3072, -49/5120, 5/7168]
+        V[1, end-1:end] = [-5/7168, 49/5120]
+
+    end
+
+    G = V[1, idx]
+    G = G ./ dx;
+
+end
+
+
+"""
+    gradNonUniform(k, ticks)
+
+Returns a m + 1 by m + 2 one-dimensional non-uniform mimetic gradient operator
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `ticks::AbstractVector`: Centers' ticks e.g. (0 0.5 1 3 5 7 9 9.5 10) (including the boundaries!)
+"""
+function gradNonUniform(k::Int, ticks::AbstractVector)
+
+    G = grad(k, length(ticks) - 2, 1)
+
+    if size(ticks, 1) == 1
+        J = diagm((G*ticks').^-1)
+    else
+        J = diagm((G*ticks).^-1)
+    end
+
+    G = J * G
+
+end
+
+#-----------------------
+# 2-D Gradient Operators
+#-----------------------
+
+"""
+    grad(k, m, dx, n, dy; dc, nc)
+
+Returns a two-dimensional mimetic gradient operator. Default is non periodic.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells in x-direction
+- `dx::T`: Step size in x-direction
+- `n::Int`: Number of cells in y-direction
+- `dy::T`: Step size in y-direction
+- `dc::NTuple{4,T}`: Dirichlet coefficients of the left, right, bottom, and top boundaries (optional)
+- `nc::NTuple{4,T}`: Neumann coefficients of the left, right, bottom, and top boundaries (optional)
+"""
+function grad(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4,T} = (1.0, 1.0, 1.0, 1.0), nc::NTuple{4,T} = (1.0, 1.0, 1.0, 1.0)) where {T}
+
+    hasbclr = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T))
+    hasbcbt = (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T))
+
+    if hasbclr
+        Gx = gradNonPeriodic(k, m, dx)
+        Im = Matrix(I, m + 2, m + 2)
+        Im = Im[2:end-1, :]
+    else
+        Gx = gradPeriodic(k, m, dx)
+        Im = Matrix(I, m, m)
+    end
+
+    if hasbcbt
+        Gy = gradNonPeriodic(k, n, dy)
+        In = Matrix(I, n + 2, n + 2)
+        In = In[2:end-1, :]
+    else
+        Gy = gradPeriodic(k, n, dy)
+        In = Matrix(I, n, n)
+    end
+
+    Sx = kron(In, Gx)
+    Sy = kron(Gy, Im)
+
+    G = [Sx; Sy];
+
+end
+
+
+"""
+    grad(k, xticks, yticks)
+
+Returns a two-dimensional non-uniform mimetic gradient operator
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `xticks::AbstractVector`: Centers' ticks (x-axis) (includint the boundaries!)
+- `yticks::AbstractVector`: Centers' ticks (y-axis) (includint the boundaries!)
+"""
+function grad(k::Int, xticks::AbstractVector, yticks::AbstractVector)
+    return gradNonUniform(k, xticks, yticks)
+end
+
+
+"""
+    gradNonPeriodic(k, m, dx, n, dy)
+
+Returns a two-dimensional non periodic mimetic gradient operator.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells in x-direction
+- `dx`: Step size in x-direction
+- `n::Int`: Number of cells in y-direction
+- `dy`: Step size in y-direction
+"""
+function gradNonPeriodic(k::Int, m::Int, dx, n::Int, dy)
+
+    Gx = gradNonPeriodic(k, m, dx)
+    Gy = gradNonPeriodic(k, n, dy)
+
+    Im = Matrix(I, m + 2, m + 2)
+    In = Matrix(I, n + 2, n + 2)
+
+    Im = Im[2:end-1, :]
+    In = In[2:end-1, :]
+
+    Sx = kron(In, Gx)
+    Sy = kron(Gy, Im)
+
+    G = [Sx; Sy];
+    
+end
+
+
+"""
+    gradPeriodic(k, m, dx, n, dy)
+
+Returns a two-dimensional periodic mimetic gradient operator.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells in x-direction
+- `dx`: Step size in x-direction
+- `n::Int`: Number of cells in y-direction
+- `dy`: Step size in y-direction
+"""
+function gradPeriodic(k::Int, m::Int, dx, n::Int, dy)
+
+    Gx = gradPeriodic(k, m, dx)
+    Gy = gradPeriodic(k, n, dy)
+
+    Im = Matrix(I, m, m)
+    In = Matrix(I, n, n)
+
+    Sx = kron(In, Gx)
+    Sy = kron(Gy, Im)
+
+    G = [Sx; Sy];
+
+end
+
+
+"""
+    gradNonUniform(k, xticks, yticks)
+
+Returns a two-dimensional non-uniform mimetic gradient operator
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `xticks::AbstractVector`: Centers' ticks (x-axis) (includint the boundaries!)
+- `yticks::AbstractVector`: Centers' ticks (y-axis) (includint the boundaries!)
+"""
+function gradNonUniform(k::Int, xticks::AbstractVector, yticks::AbstractVector)
+
+    Gx = gradNonUniform(k, xticks)
+    Gy = gradNonUniform(k, yticks)
+
+    m = size(Gx, 2)
+    n = size(Gy, 2)
+
+    Im = Matrix(I, m, m)
+    In = Matrix(I, n, n)
+
+    Im = Im[2:end-1, :]
+    In = In[2:end-1, :]
+
+    Sx = kron(In, Gx)
+    Sy = kron(Gy, Im)
+
+    G = [Sx; Sy]
+
 end
\ No newline at end of file
diff --git a/julia/MOLE.jl/src/Operators/laplacian.jl b/julia/MOLE.jl/src/Operators/laplacian.jl
index 3153a3bb..6aed4c08 100644
--- a/julia/MOLE.jl/src/Operators/laplacian.jl
+++ b/julia/MOLE.jl/src/Operators/laplacian.jl
@@ -4,19 +4,53 @@
     See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details.
 =#
 
+# -----------------------
+# 1-D Laplacian Operators
+# -----------------------
+
 """
-    lap(k, m, dx)
+    lap(k, m, dx; dc, nc)
 
-Returns a m+2 by m+2 one-dimensional mimetic laplacian operator.
+Returns a m+2 by m+2 one-dimensional mimetic laplacian operator. Default is non periodic.
 
 # Arguments
 - `k::Int`: Order of accuracy
 - `m::Int`: Number of cells
-- `dx`: Step size
+- `dx::T`: Step size
+- `dc::NTuple{2,T}`: Dirichlet coefficients of the left and right boundaries (optional)
+- `nc::NTuple{2,T}`: Neumann coefficients of the left and right boundaries (optional)
 """
-function lap(k::Int, m::Int, dx)
-    D = div(k,m,dx)
-    G = grad(k,m,dx)
+function lap(k::Int, m::Int, dx; dc::NTuple{2,T} = (1.0, 1.0), nc::NTuple{2,T} = (1.0, 1.0)) where {T}
+    D = div(k, m, dx, dc=dc, nc=nc)
+    G = grad(k, m, dx, dc=dc, nc=nc)
 
     L = D*G;
+end
+
+
+# -----------------------
+# 2-D Laplacian Operators
+# -----------------------
+
+"""
+    lap(k, m, dx, n, dy; dc, nc)
+
+Returns a two-dimensional mimetic laplacian operator. Default is non periodic.
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells in x-direction
+- `dx::T`: Step size in x-direction
+- `n::Int`: Number of cells in y-direction
+- `dy::T`: Step size in y-direction
+- `dc::NTuple{4,T}`: Dirichlet coefficients of the left and right boundaries (optional)
+- `nc::NTuple{4,T}`: Neumann coefficients of the left and right boundaries (optional)
+"""
+function lap(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4}=(1.0, 1.0, 1.0, 1.0), nc::NTuple{4}=(1.0, 1.0, 1.0, 1.0))
+
+    D = div(k, m, dx, n, dy, dc=dc, nc=nc)
+    G = grad(k, m, dx, n, dy, dc=dc, nc=nc)
+
+    L = D*G
+
 end
\ No newline at end of file
diff --git a/julia/MOLE.jl/test/BCs/scalarBC.jl b/julia/MOLE.jl/test/BCs/scalarBC.jl
index 1f2f5160..179a0672 100644
--- a/julia/MOLE.jl/test/BCs/scalarBC.jl
+++ b/julia/MOLE.jl/test/BCs/scalarBC.jl
@@ -1,3 +1,5 @@
+using Plots
+
 @testset "addScalarBC! (1D) tests" begin
 
     @testset "no BCs test" begin
@@ -142,4 +144,225 @@
         @test norm(b2 - b_ref) ≤ 1e-12
     end
 
+end
+
+@testset "addScalarBC! (2D) tests" begin
+    @testset "Periodic BCs (no BCs) test" begin
+        # Problem description
+        k = 2
+        m = 5
+        n = m + 1
+        dx = 1 / m
+        dy = 1 / n
+
+        # Boundary conditions
+        dc = (0.0, 0.0, 0.0, 0.0)
+        nc = (0.0, 0.0, 0.0, 0.0)
+        v = ([0.0], [0.0], [0.0], [0.0])
+        bc = BCs.ScalarBC2D(dc, nc, v)
+
+        # Non trivial sparse A and b
+        A = sparse(rand(m * n, m * n))
+        b = rand(m * n,)
+
+        A2, b2 = BCs.addScalarBC!(A, b, bc, k, m, dx, n, dy)
+
+        # Periodic BCs leave system unchanged
+        @test A2 == A
+        @test b2 == b
+    end
+
+    @testset "Dirichlet/Neumann BC with grad and dense matrix reference" begin
+        # Problem description
+        k = 2
+        m = 5
+        n = m + 1
+        numC = (m + 2) * (n + 2)
+        numF = (m + 1) * n + m * (n + 1)
+        dx = 0.1
+        dy = 0.4
+
+        # Boundary conditions
+        dc = (3.0, 2.0, 1.0, 4.0)
+        nc = (1.0, 8.0, 3.0, 6.0)
+        v = (rand(n,), rand(n,), rand(m + 2,), rand(m + 2,))
+        bc = BCs.ScalarBC2D(dc, nc, v)
+
+        # Non trivial sparse A and b
+        A_dense = rand(numC, numC)
+        A = sparse(A_dense)
+        b = rand(numC,)
+
+        A2, b2 = BCs.addScalarBC!(A, b, bc, k, m, dx, n, dy)
+
+        # Expected reference (dense build, no requirement that result is dense)
+        A_ref = copy(A_dense)
+        b_ref = copy(b)
+
+        # zero boundary rows and RHS entries
+        Im = Matrix(I, m + 2, m + 2)
+        In = Matrix(I, n + 2, n + 2)
+        bdry_x = zeros(Int64, m + 2, m + 2)
+        bdry_x[1, 1] = 1
+        bdry_x[end, end] = 1
+        bdry_y = zeros(Int64, n + 2, n + 2)
+        bdry_y[1, 1] = 1
+        bdry_y[end, end] = 1
+        bdry_rows = findall(vec(any(kron(In, bdry_x) .+ kron(bdry_y, Im) .!= 0, dims = 2)))
+
+        A_ref[bdry_rows, :] .= 0.0
+        b_ref[bdry_rows] .= 0.0
+
+        Gl = Matrix(Operators.grad(k, m, dx))
+        Gr = Matrix(Operators.grad(k, m, dx))
+        Gb = Matrix(Operators.grad(k, n, dy))
+        Gt = Matrix(Operators.grad(k, n, dy))
+
+        Al = zeros(Float64, m + 2, m + 2)
+        Ar = zeros(Float64, m + 2, m + 2)
+        Ab = zeros(Float64, n + 2, n + 2)
+        At = zeros(Float64, n + 2, n + 2)
+
+        Al[1, 1] = dc[1]
+        Ar[end, end] = dc[2]
+        Ab[1, 1] = dc[3]
+        At[end, end] = dc[4]
+
+        Inn = copy(In)
+        Inn[1, 1] = 0
+        Inn[end, end] = 0
+        
+        Bl = zeros(Float64, m + 2, m + 1)
+        Br = zeros(Float64, m + 2, m + 1)
+        Bb = zeros(Float64, n + 2, n + 1)
+        Bt = zeros(Float64, n + 2, n + 1)
+
+        Bl[1, 1] = -nc[1]
+        Br[end, end] = nc[2]
+        Bb[1, 1] = -nc[3]
+        Bt[end, end] = nc[4]
+
+        Al = Al + Bl * Gl
+        Ar = Ar + Br * Gr
+        Ab = Ab + Bb * Gb
+        At = At + Bt * Gt
+
+        Al = kron(Inn, Al)
+        Ar = kron(Inn, Ar)
+        Ab = kron(Ab, Im)
+        At = kron(At, Im)
+
+        A_ref .+= Al .+ Ar .+ Ab .+ At
+
+        bdry_v = zeros(size(bdry_rows,))
+        bdry_v[1:length(v[3])] = v[3]
+        midd = []
+        for i = 1:length(v[1]); midd = [midd; v[1][i]; v[2][i]]; end
+        bdry_v[length(v[3])+1:length([v[1];v[2];v[3]])] = midd
+        bdry_v[end-length(v[4])+1:end] = v[4]
+        b_ref[bdry_rows] = bdry_v
+
+
+        @test isapprox(norm(A2 - sparse(A_ref)), 0.0; rtol=1e-12, atol=1e-12)
+        @test isapprox(norm(b2 - b_ref), 0.0; rtol=1e-12, atol=1e-12)
+
+    end
+
+    @testset "Dirichlet/Neumann BC with sparse format reference" begin
+        # Problem description
+        k = 2
+        m = 5
+        n = m + 1
+        dx = 0.21
+        dy = 0.3
+        numC = (m + 2) * (n + 2)
+
+        # Nontrivial sparse A and b
+        A_dense = rand(numC, numC)
+        A = sparse(A_dense)
+        b = rand(numC,)
+
+        # Boundary conditions
+        dc = (1.0, 3.3, 4.0, 2.0)
+        nc = (1.0, 3.3, 4.0, 2.0)
+        v = (rand(n,), rand(n,), rand(m + 2,), rand(m + 2,))
+        bc = BCs.ScalarBC2D(dc, nc, v)
+
+        A2, b2 = BCs.addScalarBC!(A, b, bc, k, m, dx, n, dy)
+
+        @test A2 isa SparseMatrixCSC
+
+        # reference sparse matrix
+        A_ref = sparse(A)
+
+        # remove existing boundary rows
+        Im = Matrix(I, m + 2, m + 2)
+        In = Matrix(I, n + 2, n + 2)
+        bdry_x = zeros(Int64, m + 2, m + 2)
+        bdry_x[1, 1] = 1
+        bdry_x[end, end] = 1
+        bdry_y = zeros(Int64, n + 2, n + 2)
+        bdry_y[1, 1] = 1
+        bdry_y[end, end] = 1
+        bdry_rows = findall(vec(any(kron(In, bdry_x) .+ kron(bdry_y, Im) .!= 0, dims = 2)))
+        
+        A_ref[bdry_rows, :] = spzeros(size(A_ref[bdry_rows, :]))
+
+        # RHS boundaries
+        b_ref = copy(b)
+        bdry_v = zeros(size(bdry_rows,))
+        bdry_v[1:length(v[3])] = v[3]
+        midd = []
+        for i = 1:length(v[1]); midd = [midd; v[1][i]; v[2][i]]; end
+        bdry_v[length(v[3])+1:length([v[1];v[2];v[3]])] = midd
+        bdry_v[end-length(v[4])+1:end] = v[4]
+        b_ref[bdry_rows] = bdry_v
+
+        # LHS boundaries
+        Gl = sparse(Operators.grad(k, m, dx))
+        Gr = sparse(Operators.grad(k, m, dx))
+        Gb = sparse(Operators.grad(k, n, dy))
+        Gt = sparse(Operators.grad(k, n, dy))
+
+        Al = spzeros(Float64, m + 2, m + 2)
+        Ar = spzeros(Float64, m + 2, m + 2)
+        Ab = spzeros(Float64, n + 2, n + 2)
+        At = spzeros(Float64, n + 2, n + 2)
+
+        Al[1, 1] = dc[1]
+        Ar[end, end] = dc[2]
+        Ab[1, 1] = dc[3]
+        At[end, end] = dc[4]
+
+        Inn = copy(In)
+        Inn[1, 1] = 0
+        Inn[end, end] = 0
+        
+        Bl = spzeros(Float64, m + 2, m + 1)
+        Br = spzeros(Float64, m + 2, m + 1)
+        Bb = spzeros(Float64, n + 2, n + 1)
+        Bt = spzeros(Float64, n + 2, n + 1)
+
+        Bl[1, 1] = -nc[1]
+        Br[end, end] = nc[2]
+        Bb[1, 1] = -nc[3]
+        Bt[end, end] = nc[4]
+
+        Al = Al + Bl * Gl
+        Ar = Ar + Br * Gr
+        Ab = Ab + Bb * Gb
+        At = At + Bt * Gt
+
+        Al = kron(Inn, Al)
+        Ar = kron(Inn, Ar)
+        Ab = kron(Ab, Im)
+        At = kron(At, Im)
+
+        A_ref .+= Al .+ Ar .+ Ab .+ At     
+
+        # pure sparse comparisons
+        @test norm(A2 - A_ref) ≤ 1e-12
+        @test norm(b2 - b_ref) ≤ 1e-12
+
+    end
 end
\ No newline at end of file
diff --git a/julia/MOLE.jl/test/Operators/divergence.jl b/julia/MOLE.jl/test/Operators/divergence.jl
index 45c294cc..2bc7aee8 100644
--- a/julia/MOLE.jl/test/Operators/divergence.jl
+++ b/julia/MOLE.jl/test/Operators/divergence.jl
@@ -1,8 +1,77 @@
-@testset "Testing divergence for order k=$k" for k=2:2:8
-    tol = 1e-10
+tol = 1e-10
+
+@testset "Testing non periodic 1-D divergence for order k=$k" for k=2:2:8
     m = 2*k+1
     D = Operators.div(k, m, 1/m)
     field = ones(m+1,1)
     sol = D*field
     @test norm(sol) < tol
+end
+
+@testset "Testing periodic 1-D divergence for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    dx = 1.0 / (m - 1)
+    dc = (0.0, 0.0)
+    nc = (0.0, 0.0)
+    D = Operators.div(k, m, dx, dc=dc, nc=nc)
+    field = ones(m, 1)
+    sol = D * field
+    @test norm(sol) < tol
+end
+
+@testset "Testing non uniform 1-D divergence for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    ticks = sort(rand(m + 1))
+    D = Operators.div(k, ticks)
+    field = ones(m + 1, 1)
+    sol = D * field
+    @test norm(sol) < tol
+end
+
+@testset "Testing non periodic 2-D divergence for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    n = m + 1
+    dx = 1.0 / m
+    dy = 1.0 / n
+    D = Operators.div(k, m, dx, n, dy)
+    field = ones((m + 1) * n + m * (n + 1), 1)
+    sol = D * field
+    @test norm(sol) < tol
+end
+
+@testset "Testing periodic 2-D divergence for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    n = m + 1
+    dx = 1.0 / (m - 1)
+    dy = 1.0 / (n - 1)
+    dc = (0.0, 0.0, 0.0, 0.0)
+    nc = (0.0, 0.0, 0.0, 0.0)
+    D = Operators.div(k, m, dx, n, dy, dc=dc, nc=nc)
+    field = ones(2*m*n, 1)
+    sol = D * field
+    @test norm(sol) < tol
+end
+
+@testset "Testing mixed periodicity 2-D divergence for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    n = m + 1
+    dx = 1.0 / (m - 1)
+    dy = 1.0 / n
+    dc = (0.0, 0.0, 1.0, 1.0)
+    nc = (0.0, 0.0, 0.0, 0.0)
+    D = Operators.div(k, m, dx, n, dy, dc=dc, nc=nc)
+    field = ones(m*n + m*(n+1), 1)
+    sol = D * field
+    @test norm(sol) < tol
+end
+
+@testset "Testing non uniform 2-D divergence for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    n = m + 1
+    xticks = sort(rand(m + 1))
+    yticks = sort(rand(n + 1))
+    D = Operators.div(k, xticks, yticks)
+    field = ones((m + 1) * n + m * (n + 1), 1)
+    sol = D * field
+    @test norm(sol) < tol
 end
\ No newline at end of file
diff --git a/julia/MOLE.jl/test/Operators/gradient.jl b/julia/MOLE.jl/test/Operators/gradient.jl
index d548f1c7..5f52d92a 100644
--- a/julia/MOLE.jl/test/Operators/gradient.jl
+++ b/julia/MOLE.jl/test/Operators/gradient.jl
@@ -1,9 +1,77 @@
 tol = 1e-10
 
-@testset "Testing gradient for order k=$k" for k=2:2:8
+@testset "Testing non periodic 1-D gradient for order k=$k" for k=2:2:8
     m = 2*k+1
     G = Operators.grad(k, m, 1/m)
     field = ones(m+2, 1)
     sol = G*field
     @test norm(sol) < tol
-end;
\ No newline at end of file
+end
+
+@testset "Testing periodic 1-D gradient for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    dx = 1.0 / (m - 1)
+    dc = (0.0, 0.0)
+    nc = (0.0, 0.0)
+    G = Operators.grad(k, m, dx, dc=dc, nc=dc)
+    field = ones(m, 1)
+    sol = G * field
+    @test norm(sol) < tol
+end
+
+@testset "Testing non uniform 1-D gradient for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    ticks = sort(rand(m + 2))
+    G = Operators.grad(k, ticks)
+    field = ones(m + 2, 1)
+    sol = G * field
+    @test norm(sol) < tol
+end
+
+@testset "Testing non periodic 2-D gradient for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    n = m + 1
+    dx = 1.0 / m
+    dy = 1.0 / n
+    G = Operators.grad(k, m, dx, n, dy)
+    field = ones((m + 2) * (n + 2), 1)
+    sol = G * field
+    @test norm(sol) < tol
+end
+
+@testset "Testing periodic 2-D gradient for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    n = m + 1
+    dx = 1.0 / m
+    dy = 1.0 / n
+    dc = (0.0, 0.0, 0.0, 0.0)
+    nc = (0.0, 0.0, 0.0, 0.0)
+    G = Operators.grad(k, m, dx, n, dy, dc=dc, nc=nc)
+    field = ones(m * n, 1)
+    sol = G * field
+    @test norm(sol) < tol
+end
+
+@testset "Testing mixed periodicity 2-D gradient for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    n = m + 1
+    dx = 1.0 / m
+    dy = 1.0 / n
+    dc = (0.0, 0.0, 1.0, 1.0)
+    nc = (0.0, 0.0, 1.0, 1.0)
+    G = Operators.grad(k, m, dx, n, dy, dc=dc, nc=nc)
+    field = ones(m * (n + 2), 1)
+    sol = G * field
+    @test norm(sol) < tol
+end
+
+@testset "Testing non uniform 2-D gradient for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    n = m + 1
+    xticks = sort(rand(m + 2))
+    yticks = sort(rand(n + 2))
+    G = Operators.grad(k, xticks, yticks)
+    field = ones((m + 2) * (n + 2), 1)
+    sol = G * field
+    @test norm(sol) < tol
+end
\ No newline at end of file
diff --git a/julia/MOLE.jl/test/Operators/laplacian.jl b/julia/MOLE.jl/test/Operators/laplacian.jl
index 7091192d..b1d2e77f 100644
--- a/julia/MOLE.jl/test/Operators/laplacian.jl
+++ b/julia/MOLE.jl/test/Operators/laplacian.jl
@@ -1,9 +1,57 @@
 tol = 1e-10
 
-@testset "Testing laplacian for order k=$k" for k=2:2:8
+@testset "Testing non periodic 1-D laplacian for order k=$k" for k=2:2:8
     m = 2*k+1
     L = Operators.lap(k, m, 1/m)
     field = ones(m+2, 1)
     sol = L*field
     @test norm(sol) < tol
-end;
\ No newline at end of file
+end
+
+@testset "Testing periodic 1-D laplacian for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    dx = 1.0 / m
+    dc = (0.0, 0.0)
+    nc = (0.0, 0.0)
+    L = Operators.lap(k, m, dx, dc=dc, nc=nc)
+    field = ones(m, 1)
+    sol = L * field
+    @test norm(sol) < tol
+end
+
+@testset "Testing non periodic 2-D laplacian for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    n = m + 1
+    dx = 1.0 / m
+    dy = 1.0 / n
+    L = Operators.lap(k, m, dx, n, dy)
+    field = ones((m + 2) * (n + 2), 1)
+    sol = L * field
+    @test norm(sol) < tol
+end
+
+@testset "Testing periodic 2-D laplacian for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    n = m + 1
+    dx = 1.0 / m
+    dy = 1.0 / n
+    dc = (0.0, 0.0, 0.0, 0.0)
+    nc = (0.0, 0.0, 0.0, 0.0)
+    L = Operators.lap(k, m, dx, n, dy, dc=dc, nc=nc)
+    field = ones(m * n, 1)
+    sol = L * field
+    @test norm(sol) < tol
+end
+
+@testset "Testing mixed periodicity 2-D laplacian for order k=$k" for k=2:2:8
+    m = 2 * k + 1
+    n = m + 1
+    dx = 1.0 / m
+    dy = 1.0 / n
+    dc = (0.0, 0.0, 1.0, 1.0)
+    nc = (0.0, 0.0, 0.0, 0.0)
+    L = Operators.lap(k, m, dx, n, dy, dc=dc, nc=nc)
+    field = ones(m * (n + 2), 1)
+    sol = L * field
+    @test norm(sol) < tol
+end
\ No newline at end of file