diff --git a/Project.toml b/Project.toml
index f7561c700..9cd88170c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -32,7 +32,7 @@ Interpolations = "0.14, 0.15"
 Latexify = "0.15, 0.16"
 ModelingToolkit = "9.17"
 OrdinaryDiffEq = "6"
-PDEBase = "0.1.12"
+PDEBase = "0.1.17"
 PrecompileTools = "1"
 RuntimeGeneratedFunctions = "0.5.9"
 SafeTestsets = "0.0.1"
diff --git a/docs/src/tutorials/schroedinger.md b/docs/src/tutorials/schroedinger.md
index c4232c694..0e834a213 100644
--- a/docs/src/tutorials/schroedinger.md
+++ b/docs/src/tutorials/schroedinger.md
@@ -18,9 +18,9 @@ V(x) = 0.0
 
 eq = [im * Dt(ψ(t, x)) ~ Dxx(ψ(t, x)) + V(x) * ψ(t, x)] # You must enclose complex equations in a vector, even if there is only one equation
 
-ψ0 = x -> sin(2pi * x)
+ψ0 = x -> ((1 + im)/sqrt(2))*sinpi(2*x)
 
-bcs = [ψ(0, x) ~ ψ0(x),
+bcs = [ψ(0, x) => ψ0(x), # Initial condition must be marked with a => operator
     ψ(t, xmin) ~ 0,
     ψ(t, xmax) ~ 0]
 
@@ -46,4 +46,6 @@ end
 gif(anim, "schroedinger.gif", fps = 10)
 ```
 
+Note that complex initial conditions are supported, but must be marked with a `=>` operator.
+
 This represents the second from ground state of a particle in an infinite quantum well, try changing the potential `V(x)`, initial conditions and BCs, it is extremely interesting to see how the wave function evolves even for nonphysical combinations. Be sure to post interesting results on the discourse!
diff --git a/test/pde_systems/schroedinger.jl b/test/pde_systems/schroedinger.jl
index 81e53442b..d8573bc3a 100644
--- a/test/pde_systems/schroedinger.jl
+++ b/test/pde_systems/schroedinger.jl
@@ -1,55 +1,90 @@
 using MethodOfLines, OrdinaryDiffEq, DomainSets, ModelingToolkit, Test
 
-@parameters t, x
-@variables ψ(..)
+@testset "Schroedinger" begin
+    @parameters t, x
+    @variables ψ(..)
 
-Dt = Differential(t)
-Dxx = Differential(x)^2
+    Dt = Differential(t)
+    Dxx = Differential(x)^2
 
-xmin = 0
-xmax = 1
+    xmin = 0
+    xmax = 1
 
-V(x) = 0.0
+    V(x) = 0.0
 
-eq = [im * Dt(ψ(t, x)) ~ (Dxx(ψ(t, x)) + V(x) * ψ(t, x))] # You must enclose complex equations in a vector, even if there is only one equation
+    eq = [im * Dt(ψ(t, x)) ~ (Dxx(ψ(t, x)) + V(x) * ψ(t, x))] # You must enclose complex equations in a vector, even if there is only one equation
 
-ψ0 = x -> sin(2 * pi * x)
+    ψ0 = x -> ((1 + im)/sqrt(2))*sinpi(2*x)
 
-bcs = [ψ(0, x) ~ ψ0(x),
-    ψ(t, xmin) ~ 0,
-    ψ(t, xmax) ~ 0]
+    bcs = [ψ(0, x) => ψ0(x), # Initial condition must be marked with a => operator
+        ψ(t, xmin) ~ 0,
+        ψ(t, xmax) ~ 0]
 
-domains = [t ∈ Interval(0, 1), x ∈ Interval(xmin, xmax)]
+    domains = [t ∈ Interval(0, 1), x ∈ Interval(xmin, xmax)]
 
-@named sys = PDESystem(eq, bcs, domains, [t, x], [ψ(t, x)])
+    @named sys = PDESystem(eq, bcs, domains, [t, x], [ψ(t, x)])
 
-disc = MOLFiniteDifference([x => 100], t)
+    disc = MOLFiniteDifference([x => 100], t)
 
-prob = discretize(sys, disc)
+    prob = discretize(sys, disc)
 
-sol = solve(prob, TRBDF2(), saveat = 0.01)
+    sol = solve(prob, TRBDF2(), saveat = 0.01)
 
-discx = sol[x]
-disct = sol[t]
+    discx = sol[x]
+    disct = sol[t]
 
-discψ = sol[ψ(t, x)]
+    discψ = sol[ψ(t, x)]
 
-analytic(t, x) = sqrt(2) * sin(2 * pi * x) * exp(-im * 4 * pi^2 * t)
+    analytic(t, x) = sqrt(2) * sin(2 * pi * x) * exp(-im * 4 * pi^2 * t) * ((1 + im)/sqrt(2))
 
-analψ = [analytic(t, x) for t in disct, x in discx]
+    analψ = [analytic(t, x) for t in disct, x in discx]
 
-for i in 1:length(disct)
-    u = abs.(analψ[i, :]) .^ 2
-    u2 = abs.(discψ[i, :]) .^ 2
+    for i in 1:length(disct)
+        u = abs.(analψ[i, :]) .^ 2
+        u2 = abs.(discψ[i, :]) .^ 2
 
-    @test u ./ maximum(u)≈u2 ./ maximum(u2) atol=1e-3
+        @test u ./ maximum(u)≈u2 ./ maximum(u2) atol=1e-2
+    end
+    #using Plots
+    
+    # anim = @animate for i in 1:length(disct)
+    #     u = analψ[i, :]
+    #     u2 = discψ[i, :]
+    #     plot(discx, [real.(u), imag.(u)], ylim = (-1.5, 1.5), title = "t = $(disct[i])", xlabel = "x", ylabel = "ψ(t,x)", label = ["re(ψ)" "im(ψ)"], legend = :topleft)
+    # end
+    # gif(anim, "schroedinger.gif", fps = 10)
 end
 
-#using Plots
-
-# anim = @animate for i in 1:length(disct)
-#     u = analψ[i, :]
-#     u2 = discψ[i, :]
-#     plot(discx, [real.(u), imag.(u)], ylim = (-1.5, 1.5), title = "t = $(disct[i])", xlabel = "x", ylabel = "ψ(t,x)", label = ["re(ψ)" "im(ψ)"], legend = :topleft)
-# end
-# gif(anim, "schroedinger.gif", fps = 10)
+@testset "Schroedinger with complex bcs" begin
+    @parameters t, x
+    @variables ψ(..)
+    
+    Dt = Differential(t)
+    Dxx = Differential(x)^2
+    
+    xmin = 0
+    xmax = 1
+    ϵ = 1e-2
+    
+    V(x) = 0.0
+    
+    eq = [(im*ϵ)*Dt(ψ(t,x)) ~ (-0.5*ϵ^2)Dxx(ψ(t,x)) + V(x)*ψ(t,x)] # You must enclose complex equations in a vector, even if there is only one equation
+    
+    ψ0 = x -> exp((im/ϵ)*1e-1*sum(x))
+    
+    bcs = [ψ(0,x) => ψ0(x),
+        ψ(t,xmin) ~ exp((im/ϵ)*(1e-1*sum(xmin) - 0.5e-2*t)),
+        ψ(t,xmax) ~ exp((im/ϵ)*(1e-1*sum(xmax) - 0.5e-2*t))]
+    
+    domains = [t ∈ Interval(0, 1), x ∈ Interval(xmin, xmax)]
+    
+    @named sys = PDESystem(eq, bcs, domains, [t, x], [ψ(t,x)])
+    
+    disc = MOLFiniteDifference([x => 300], t)
+    
+    prob = discretize(sys, disc)
+    
+    sol = solve(prob, TRBDF2(), saveat = 0.01)
+
+    @test SciMLBase.successful_retcode(sol)
+end 
\ No newline at end of file