diff --git a/Project.toml b/Project.toml
index de5b02fbf..b765faf98 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,7 +13,6 @@ LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
diff --git a/docs/src/examples.md b/docs/src/examples.md
index d3fe4d5ec..237e99a52 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -2,8 +2,12 @@
 
 ## Heat equation
 
+### Dirichlet boundary conditions
+
 ```julia
 using ModelingToolkit, DiffEqOperators
+# Method of Manufactured Solutions: exact solution
+u_exact = (x,t) -> exp.(-t) * cos.(x)
 
 # Parameters, variables, and derivatives
 @parameters t x
@@ -13,13 +17,13 @@ using ModelingToolkit, DiffEqOperators
 
 # 1D PDE and boundary conditions
 eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
-bcs = [u(0,x) ~ -x*(x-1)*sin(x),
-        u(t,0) ~ 0.0,
-        u(t,1) ~ 0.0]
+bcs = [u(0,x) ~ cos(x),
+        u(t,0) ~ exp(-t),
+        u(t,1) ~ exp(-t) * cos(1)]
 
 # Space and time domains
 domains = [t ∈ IntervalDomain(0.0,1.0),
-            x ∈ IntervalDomain(0.0,1.0)]
+        x ∈ IntervalDomain(0.0,1.0)]
 
 # PDE system
 pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
@@ -34,5 +38,122 @@ prob = discretize(pdesys,discretization)
 
 # Solve ODE problem
 using OrdinaryDiffEq
-sol = solve(prob,Tsit5(),saveat=0.1)
+sol = solve(prob,Tsit5(),saveat=0.2)
+
+# Plot results and compare with exact solution
+x = prob.space[2]
+t = sol.t
+
+using Plots
+plt = plot()
+for i in 1:length(t)
+    plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]),label="Numerical, t=$(t[i])")
+    scatter!(x, u_exact(x, t[i]),label="Exact, t=$(t[i])")
+end
+display(plt)
+```
+### Neumann boundary conditions
+
+```julia
+using ModelingToolkit, DiffEqOperators
+# Method of Manufactured Solutions: exact solution
+u_exact = (x,t) -> exp.(-t) * cos.(x)
+
+# Parameters, variables, and derivatives
+@parameters t x
+@variables u(..)
+@derivatives Dt'~t
+@derivatives Dx'~x
+@derivatives Dxx''~x
+
+# 1D PDE and boundary conditions
+eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
+bcs = [u(0,x) ~ cos(x),
+        Dx(u(t,0)) ~ 0.0,
+        Dx(u(t,1)) ~ -exp(-t) * sin(1)]
+
+# Space and time domains
+domains = [t ∈ IntervalDomain(0.0,1.0),
+        x ∈ IntervalDomain(0.0,1.0)]
+
+# PDE system
+pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+
+# Method of lines discretization
+# Need a small dx here for accuracy
+dx = 0.01
+order = 2
+discretization = MOLFiniteDifference(dx,order)
+
+# Convert the PDE problem into an ODE problem
+prob = discretize(pdesys,discretization)
+
+# Solve ODE problem
+using OrdinaryDiffEq
+sol = solve(prob,Tsit5(),saveat=0.2)
+
+# Plot results and compare with exact solution
+x = prob.space[2]
+t = sol.t
+
+using Plots
+plt = plot()
+for i in 1:length(t)
+    plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]),label="Numerical, t=$(t[i])")
+    scatter!(x, u_exact(x, t[i]),label="Exact, t=$(t[i])")
+end
+display(plt)
+```
+
+### Robin boundary conditions
+
+```julia
+using ModelingToolkit, DiffEqOperators 
+# Method of Manufactured Solutions
+u_exact = (x,t) -> exp.(-t) * sin.(x)
+
+# Parameters, variables, and derivatives
+@parameters t x
+@variables u(..)
+@derivatives Dt'~t
+@derivatives Dx'~x
+@derivatives Dxx''~x
+
+# 1D PDE and boundary conditions
+eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
+bcs = [u(0,x) ~ sin(x),
+        u(t,-1.0) + 3Dx(u(t,-1.0)) ~ exp(-t) * (sin(-1.0) + 3cos(-1.0)),
+        u(t,1.0) + Dx(u(t,1.0)) ~ exp(-t) * (sin(1.0) + cos(1.0))]
+
+# Space and time domains
+domains = [t ∈ IntervalDomain(0.0,1.0),
+        x ∈ IntervalDomain(-1.0,1.0)]
+
+# PDE system
+pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+
+# Method of lines discretization
+# Need a small dx here for accuracy
+dx = 0.05
+order = 2
+discretization = MOLFiniteDifference(dx,order)
+
+# Convert the PDE problem into an ODE problem
+prob = discretize(pdesys,discretization)
+
+# Solve ODE problem
+using OrdinaryDiffEq
+sol = solve(prob,Tsit5(),saveat=0.2)
+
+# Plot results and compare with exact solution
+x = prob.space[2]
+t = sol.t
+
+using Plots
+plt = plot()
+for i in 1:length(t)
+    plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]),label="Numerical, t=$(t[i])")
+    scatter!(x, u_exact(x, t[i]),label="Exact, t=$(t[i])")
+end
+display(plt)
 ```
\ No newline at end of file
diff --git a/src/MOL_discretization.jl b/src/MOL_discretization.jl
index 50cd63e40..32a7acf33 100644
--- a/src/MOL_discretization.jl
+++ b/src/MOL_discretization.jl
@@ -94,7 +94,7 @@ function get_bcs(bcs,tdomain,domain)
         else
             throw(BoundaryConditionError(
                 "Boundary condition not recognized. "
-                * "Should be applied at t=0, x=$(domain.lower), or x=$(domain.upper)"
+                * "Should be applied at t=$(tdomain.lower), x=$(domain.lower), or x=$(domain.upper)"
             ))
         end
         # Create value