diff --git a/src/MOLFiniteDifference/MOL_discretization.jl b/src/MOLFiniteDifference/MOL_discretization.jl
index c12a0059c..054148c36 100644
--- a/src/MOLFiniteDifference/MOL_discretization.jl
+++ b/src/MOLFiniteDifference/MOL_discretization.jl
@@ -140,7 +140,7 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
 
             bclocs = map(e->substitute.(indvars,e),edgevals) # location of the boundary conditions e.g. (t,0.0,y)
             edgemaps = Dict(bclocs .=> [spacevals[e...] for e in edges])
-            initmaps = depvars
+            
             if t != nothing
                 initmaps = substitute.(depvars,[t=>tspan[1]])
             end
@@ -187,13 +187,31 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
             # Generate initial conditions and bc equations
             for bc in pdesys.bcs
                 bcdepvar = first(get_depvars(bc.lhs, depvar_ops))
+                # Make sure the bcdepvar is one of the depvars
                 if any(u->isequal(operation(u),operation(bcdepvar)),depvars)
+                    # Check whether the bc is an initial condition
+                    is_ic = false
                     if t != nothing && operation(bc.lhs) isa Sym && !any(x -> isequal(x, t.val), arguments(bc.lhs))
+                        # Look for u(t,0) case
+                        depvar = bc.lhs # e.g. u(t,0)
+                        bcop = identity
+                        is_ic = true
+                    elseif t != nothing && operation(bc.lhs) isa Differential
+                        # Look for Dt(u(t,0)) case
+                        depvar = arguments(bc.lhs)[1] # e.g. u(t,0)
+                        if operation(depvar) isa Sym && !any(x -> isequal(x, t.val), arguments(depvar))
+                            bcop = operation(bc.lhs) # e.g. Dt
+                            is_ic = true
+                        end
+                    end
+                    if is_ic
                         # initial condition
-                        # Assume in the form `u(...) ~ ...` for now
-                        i = findfirst(isequal(bc.lhs),initmaps)
+                        # Assume in the form `u(...) ~ ...` or `Dt(u(...)) ~ ...` for now
+                        # find the index of the depvar
+                        i = findfirst(isequal(depvar),initmaps)
                         if i !== nothing
-                            push!(u0,vec(depvarsdisc[i] .=> substitute.((bc.rhs,),gridvals)))
+                            # bcop is identity for u(...) and Dt for Dt(u(...))
+                            push!(u0,vec(bcop.(depvarsdisc[i]) .=> substitute.((bc.rhs,),gridvals)))
                         end
                     else
                         # Algebraic equations for BCs
@@ -207,7 +225,7 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
                             # Replace symbol in the bc lhs with the spatial discretized term
                             lhs = substitute(lhs,depvarbcmaps[i])
                             rhs = substitute.((bc.rhs,),edgemaps[bcargs])
-                            lhs = lhs isa Vector ? lhs : [lhs] # handle 1D
+                            lhs = lhs isa Array ? lhs : [lhs] # handle 1D
                             push!(bceqs,lhs .~ rhs)
                         end
                     end
@@ -375,6 +393,8 @@ end
 
 function SciMLBase.discretize(pdesys::ModelingToolkit.PDESystem,discretization::DiffEqOperators.MOLFiniteDifference)
     sys, tspan = SciMLBase.symbolic_discretize(pdesys,discretization)
+    # Lower order for higher-order ODEs (e.g. wave equation)
+    sys = sys isa NonlinearSystem ? sys : ode_order_lowering(sys)
     if tspan == nothing
         return prob = NonlinearProblem(sys, ones(length(sys.states)))
     else
diff --git a/test/MOL/MOL_1D_Wave.jl b/test/MOL/MOL_1D_Wave.jl
new file mode 100644
index 000000000..c0a5c86f6
--- /dev/null
+++ b/test/MOL/MOL_1D_Wave.jl
@@ -0,0 +1,55 @@
+# 1D wave equation problems
+
+# Packages and inclusions
+using ModelingToolkit,DiffEqOperators,LinearAlgebra,Test,OrdinaryDiffEq, DomainSets
+using ModelingToolkit: Differential
+
+# Tests
+# @testset "Test 00: Dtt(u(t,x)) ~ Dxx(u(t,x))" begin
+    # Method of Manufactured Solutions
+    u_exact = (x,t) -> exp.(-t) * cos.(x)
+
+    # Parameters, variables, and derivatives
+    @parameters t x
+    @variables u(..)
+    Dt = Differential(t)
+    Dtt = Differential(t)^2
+    Dxx = Differential(x)^2
+
+    # 1D PDE and boundary conditions
+    eq  = Dtt(u(t,x)) ~ -Dxx(u(t,x))
+    bcs = [u(0,x) ~ cos(x),
+           Dt(u(0,x)) ~ -cos(x),
+           u(t,0) ~ exp(-t),
+           u(t,Float64(π)) ~ -exp(-t)]
+
+    # Space and time domains
+    domains = [t ∈ Interval(0.0,1.0),
+               x ∈ Interval(0.0,Float64(π))]
+
+    # PDE system
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])
+
+    # Method of lines discretization
+    dx = range(0.0,Float64(π),length=30)
+    order = 2
+    discretization = MOLFiniteDifference([x=>dx],t)
+
+    # Convert the PDE problem into an ODE problem
+    prob = discretize(pdesys,discretization)
+    sys, tspan = SciMLBase.symbolic_discretize(pdesys,discretization)
+    # Solve ODE problem
+    sol = solve(prob,Tsit5(),saveat=0.1)
+    x_sol = dx[2:end-1]
+    t_sol = sol.t
+
+    # Test against exact solution
+    # for i in 1:length(sol)
+    i = 10
+        exact = u_exact(x_sol, t_sol[i])
+        # non-differential states only (ode_sys_lowering creates additional differential states)
+        u_approx = sol.u[i][1:28]
+        @test all(isapprox.(u_approx, exact, atol=0.01))
+    # end
+# end
+
diff --git a/test/MOL/MOL_NonlinearProblem.jl b/test/MOL/MOL_NonlinearProblem.jl
index 88efcfe57..1044f8afe 100644
--- a/test/MOL/MOL_NonlinearProblem.jl
+++ b/test/MOL/MOL_NonlinearProblem.jl
@@ -1,4 +1,4 @@
-# 1D diffusion problem
+# 1D nonlinear problem
 
 # Packages and inclusions
 using ModelingToolkit,DiffEqOperators,LinearAlgebra,Test
diff --git a/test/runtests.jl b/test/runtests.jl
index e58ca672c..2e9c47e83 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -48,6 +48,7 @@ if GROUP == "All" || GROUP == "MOLFiniteDifference"
     @time @safetestset "MOLFiniteDifference Interface: 1D HigherOrder" begin include("MOL/MOL_1D_HigherOrder.jl") end
     @time @safetestset "MOLFiniteDifference Interface: 1D Partial DAE" begin include("MOL/MOL_1D_PDAE.jl") end
     @time @safetestset "MOLFiniteDifference Interface: Stationary Nonlinear Problems" begin include("MOL/MOL_NonlinearProblem.jl") end
+    @time @safetestset "MOLFiniteDifference Interface: 1D Wave Equation" begin include("MOL/MOL_1D_Wave.jl") end
 end
 
 if GROUP == "All" || GROUP == "Misc"