diff --git a/julia/MOLE.jl/examples/hyperbolic1D.jl b/julia/MOLE.jl/examples/hyperbolic1D.jl
new file mode 100644
index 00000000..0c60354e
--- /dev/null
+++ b/julia/MOLE.jl/examples/hyperbolic1D.jl
@@ -0,0 +1,71 @@
+#=
+Julia implementation of the 1D Advection equation with periodic boundary conditions
+from MATLAB MOLE example
+=#
+
+using Printf
+using Plots
+import MOLE: Operators, BCs
+
+# Domain limits
+west = 0;
+east = 1;
+
+a  = 1; #velocity
+
+
+k  = 2; #Operator's order of accuracy
+m  = 50; # number of cells
+dx = (east - west) / m;
+
+t  = 1; #Simulation time
+dt = dx/abs(a); #CFL condition for explicit schemes
+
+D  = Operators.div(k,m,dx); #1D Mimetic divergence operator
+I  = Operators.interpol(m, 0.5); #1D 2nd order interpolator
+
+# 1D Staggered grid
+grid = [west, west+dx/2:dx:east-dx/2, east];
+
+# Initial Condition
+q = 2 * π .* grid
+println(typeof(q))
+U = sin.(q)';
+
+#Periodic Boundary Condition imposed on the divergence operator
+D[1, 2]       = 1/(2*dx);
+D[1, end-1]   = -1/(2*dx);
+D[end, 1]     = -1/(2*dx);
+D[end, end-1] = -1/(2*dx);
+
+#Premultiply out of the time loop (since it does not change)
+D = -a*dt*2*D*I;
+
+#=
+        We could have also defined 
+        D = -a*dt*2*I*D
+        if the grid was defined as:
+        grid = west : dx :east (nodal)
+=#
+
+U2 = U + D/2*U #Compute one step using Euler's method
+
+
+#Time integration loop
+for i in 1: t/dt
+    #Plot approximation
+    plot(grid, U2, label="Approximated")
+    #Plot exact solution
+    plot!(grid, sin.(2*π .*(grid - a*i*dt)), label="exact")
+    #Plot attributes
+    title!("1D Advection Equation with Periodic BC");
+    xlabel!("x");
+    ylabel!("u(x,t)");
+    axis!([west east -1.5 1.5]);
+    sleep(0.04);
+    U3 = U + D*U2; #Compute next step using leapfrog scheme
+    U  = U2;
+    U2 = U3;
+end
+
+
diff --git a/julia/MOLE.jl/src/Operators/Operators.jl b/julia/MOLE.jl/src/Operators/Operators.jl
index 20491d46..7af4cfb7 100644
--- a/julia/MOLE.jl/src/Operators/Operators.jl
+++ b/julia/MOLE.jl/src/Operators/Operators.jl
@@ -3,5 +3,6 @@ module Operators
 include("gradient.jl")
 include("divergence.jl")
 include("laplacian.jl")
+include("interpol.jl")
 
-end # module Operators
\ No newline at end of file
+end # module Operators
diff --git a/julia/MOLE.jl/src/Operators/interpol.jl b/julia/MOLE.jl/src/Operators/interpol.jl
new file mode 100644
index 00000000..68c9d401
--- /dev/null
+++ b/julia/MOLE.jl/src/Operators/interpol.jl
@@ -0,0 +1,40 @@
+using SparseArrays
+#=
+        SPDX-License-Identifier: GPL-3.0-or-later
+        © 2008-2024 San Diego State University Research Foundation (SDSURF).
+        See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details.
+=#
+
+function interpol(m::Int,c::Float64)
+    """
+    implementation of MATLAB interpol.m operator for Julia
+    INPUTS: 
+       'm::Int'   : number of cells
+       'c::Float' : left interpolation coefficient
+    OUTPUTS:
+       I : a (m+1)×(m+2) one-dimensional interpolator of 2nd-order
+    """
+
+    #Assertions:
+    @assert m>=4 ["m >= 4"];
+    @assert c >= 0 && c <= 1 ["0 <= c <= 1"];
+
+    #Dimensions of I:
+    n_rows = m+1;
+    n_cols = m+2;
+
+    I = spzeros(n_rows, n_cols);
+
+    I[1,1] = 1;
+    I[end,end]=1;
+
+    #Average between two continuous cells
+    avg = [c 1-c];
+
+    j = 2;
+    for i = 2: n_rows - 1
+        I[i, j:j+2-1] = avg;
+        j = j + 1;
+    end
+end
+