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/DiffEqOperators.jl b/src/DiffEqOperators.jl
index fe9918484..0af00efb2 100644
--- a/src/DiffEqOperators.jl
+++ b/src/DiffEqOperators.jl
@@ -21,6 +21,9 @@ include("jacvec_operators.jl")
 ### Utilities
 include("utils.jl")
 
+### Exceptions
+include("exceptions.jl")
+
 ### Boundary Padded Arrays
 include("boundary_padded_arrays.jl")
 
@@ -63,4 +66,5 @@ export compose, decompose, perpsize
 
 export GhostDerivativeOperator
 export MOLFiniteDifference
+export BoundaryConditionError
 end # module
diff --git a/src/MOL_discretization.jl b/src/MOL_discretization.jl
index fcdf49508..2c2552ec5 100644
--- a/src/MOL_discretization.jl
+++ b/src/MOL_discretization.jl
@@ -6,29 +6,129 @@ struct MOLFiniteDifference{T} <: DiffEqBase.AbstractDiscretization
     MOLFiniteDifference(args...;order=2) = new{typeof(args[1])}(args[1],order)
 end
 
+function throw_bc_err(bc)
+    throw(BoundaryConditionError(
+        "Could not read boundary condition '$(bc.lhs) ~ $(bc.rhs)'"
+    ))
+end
 # Get boundary conditions from an array
+# The returned boundary conditions will be the three coefficients (α, β, γ)
+# for a RobinBC of the form α*u(0) + β*u'(0) = γ
+# These are created automatically from the form of the boundary condition that is
+# passed, i.e.
+# u(t, 0) ~ γ                           ---> (1, 0, γ)
+# Dx(u(t, 0)) ~ γ                       ---> (0, 1, γ)
+# α * u(t, 0)) + β * Dx(u(t, 0)) ~ γ    ---> (α, β, γ)
+# In general, all of α, β, and γ could be Expr (i.e. functions of t)
 function get_bcs(bcs,tdomain,domain)
     lhs_deriv_depvars_bcs = Dict()
     num_bcs = size(bcs,1)
     for i = 1:num_bcs
-        var = operation(bcs[i].lhs)
-        if var isa Sym
-            var = var.name
-            if !haskey(lhs_deriv_depvars_bcs,var)
-                lhs_deriv_depvars_bcs[var] = Array{Expr}(undef,3)
-            end
-            j = 0
-            if isequal(bcs[i].lhs.args[1],tdomain.lower) # u(t=0,x)
-                j = 1
-            elseif isequal(bcs[i].lhs.args[2],domain.lower) # u(t,x=x_init)
-                j = 2
-            elseif isequal(bcs[i].lhs.args[2],domain.upper) # u(t,x=x_final)
-                j = 3
+        lhs = bcs[i].lhs
+        # Extract the variable from the lhs
+        if operation(lhs) isa Sym
+            # Dirichlet boundary condition
+            var = nameof(operation(lhs))
+            α = 1.0
+            β = 0.0
+            bc_args = lhs.args
+        elseif operation(lhs) isa Differential
+            # Neumann boundary condition
+            # Check that we don't have a second-order derivative in the
+            # boundary condition, by checking that the argument is a Sym
+            @assert operation(lhs.args[1]) isa Sym throw_bc_err(bcs[i])
+            var = nameof(operation(lhs.args[1]))
+            α = 0.0
+            β = 1.0
+            bc_args = lhs.args[1].args
+        elseif operation(lhs) isa typeof(+)
+            # Robin boundary condition
+            lhs_l, lhs_r = lhs.args
+            # Left side of the expression should be Sym or α * Sym
+            if operation(lhs_l) isa Sym
+                α = 1.0
+                var_l = nameof(operation(lhs_l))
+                bc_args_l = lhs_l.args
+            elseif operation(lhs_l) isa typeof(*)
+                α = lhs_l.args[1]
+                # Convert α to a Float64 if it is an Int, leave unchanged otherwise
+                α = α isa Int ? Float64(α) : α
+                @assert operation(lhs_l.args[2]) isa Sym throw_bc_err(bcs[i])
+                var_l = nameof(operation(lhs_l.args[2]))
+                bc_args_l = lhs_l.args[2].args
+            else
+                throw_bc_err(bcs[i])
             end
-            if !(bcs[i].rhs isa ModelingToolkit.Symbolic)
-                lhs_deriv_depvars_bcs[var][j] = :(var=$(bcs[i].rhs))
+            # Right side of the expression should be Differential or β * Differential
+            if operation(lhs_r) isa Differential
+                # Check that we don't have a second-order derivative in the
+                # boundary condition
+                @assert operation(lhs_r.args[1]) isa Sym throw_bc_err(bcs[i])
+                β = 1.0
+                var_r = nameof(operation(lhs_r.args[1]))
+                bc_args_r = lhs_r.args[1].args
+            elseif operation(lhs_r) isa typeof(*)
+                β = lhs_r.args[1]
+                # Convert β to a Float64 if it is an Int, leave unchanged otherwise
+                β = β isa Int ? Float64(β) : β
+                # Check that the bc is a derivative
+                @assert operation(lhs_r.args[2]) isa Differential throw_bc_err(bcs[i])
+                # But not second order (argument should be a Sym)
+                @assert operation(lhs_r.args[2].args[1]) isa Sym throw_bc_err(bcs[i])
+                var_r = nameof(operation(lhs_r.args[2].args[1]))
+                bc_args_r = lhs_r.args[2].args[1].args
             else
-                lhs_deriv_depvars_bcs[var][j] = toexpr(bcs[i].rhs)
+                throw_bc_err(bcs[i])
+            end
+            # Check var and bc_args are the same in lhs and rhs, and if so assign 
+            # the unique value
+            @assert var_l == var_r throw(BoundaryConditionError(
+                "mismatched variables '$var_l' and '$var_r' "
+                * "in Robin BC '$(bcs[i].lhs) ~ $(bcs[i].rhs)'"
+            ))
+            var = var_l
+            @assert bc_args_l == bc_args_r throw(BoundaryConditionError(
+                "mismatched args $bc_args_l and $bc_args_r "
+                * "in Robin BC '$(bcs[i].lhs) ~ $(bcs[i].rhs)'"
+            ))
+            bc_args = bc_args_l
+        else
+            throw_bc_err(bcs[i])
+        end
+        if !haskey(lhs_deriv_depvars_bcs,var)
+            # Initialize dict of boundary conditions for this variable
+            lhs_deriv_depvars_bcs[var] = Dict()
+        end
+        # Create key
+        if isequal(bc_args[1],tdomain.lower) # u(t=0,x)
+            key = "ic"
+        elseif isequal(bc_args[2],domain.lower) # u(t,x=x_init)
+            key = "left bc"
+        elseif isequal(bc_args[2],domain.upper) # u(t,x=x_final)
+            key = "right bc"
+        else
+            throw(BoundaryConditionError(
+                "Boundary condition '$(bcs[i].lhs) ~ $(bcs[i].rhs)' could not be read. "
+                * "BCs should be applied at t=$(tdomain.lower), "
+                * "x=$(domain.lower), or x=$(domain.upper)"
+            ))
+        end
+        # Create value
+        γ = toexpr(bcs[i].rhs)
+        # Assign
+        if key == "ic"
+            # Initial conditions always take the form u(0, x) ~ γ
+            lhs_deriv_depvars_bcs[var][key] = γ
+        else
+            # Boundary conditions can be more general
+            lhs_deriv_depvars_bcs[var][key] = (toexpr(α), toexpr(β), γ)
+        end
+    end
+    # Check each variable got all boundary conditions
+    for var in keys(lhs_deriv_depvars_bcs)
+        for key in ["ic", "left bc", "right bc"]
+            if !haskey(lhs_deriv_depvars_bcs[var], key)
+                throw(BoundaryConditionError("missing $key for $var"))
             end
         end
     end
@@ -43,7 +143,7 @@ end
 #       E.g., Dx(u(t,x))=v(t,x)*Dx(u(t,x)), v(t,x)=t*x
 #            =>  Dx(u(t,x))=t*x*Dx(u(t,x))
 
-function discretize_2(input,deriv_order,approx_order,dx,X,len,
+function discretize_2(input,deriv_order,approx_order,dx,X,len_of_indep_vars,
                       deriv_var,dep_var_idx,indep_var_idx)
     if !(input isa ModelingToolkit.Symbolic)
         return :($(input))
@@ -54,7 +154,7 @@ function discretize_2(input,deriv_order,approx_order,dx,X,len,
             if haskey(indep_var_idx,var) # ind. var.
                 if var != :(t)
                     i = indep_var_idx[var]
-                    expr = :($X[$i][2:$len[$i]-1])
+                    expr = :($X[$i][2:$len_of_indep_vars[$i]-1])
                 else
                     expr = :(t)
                 end
@@ -68,11 +168,11 @@ function discretize_2(input,deriv_order,approx_order,dx,X,len,
                     # TODO: approx_order and forward/backward should be
                     #       input parameters of each derivative
                     approx_order = 1
-                    L = UpwindDifference(deriv_order,approx_order,dx[i],len[i]-2,-1)
-                    expr = :(-1*($L*all_bcs[$j]*u[:,$j]))
+                    L = UpwindDifference(deriv_order,approx_order,dx[i],len_of_indep_vars[i]-2,-1)
+                    expr = :(-1*($L*Q[$j]*u[:,$j]))
                 elseif deriv_order == 2
-                    L = CenteredDifference(deriv_order,approx_order,dx[i],len[i]-2)
-                    expr = :($L*all_bcs[$j]*u[:,$j])
+                    L = CenteredDifference(deriv_order,approx_order,dx[i],len_of_indep_vars[i]-2)
+                    expr = :($L*Q[$j]*u[:,$j])
                 end
             end
             return expr
@@ -80,18 +180,18 @@ function discretize_2(input,deriv_order,approx_order,dx,X,len,
             var = nameof(input.op.x)
             push!(deriv_var,var)
             return discretize_2(input.args[1],deriv_order+1,approx_order,dx,X,
-                                len,deriv_var,dep_var_idx,indep_var_idx)
+                                len_of_indep_vars,deriv_var,dep_var_idx,indep_var_idx)
         else
             name = nameof(operation(input))
             if size(input.args,1) == 1
                 aux = discretize_2(input.args[1],deriv_order,approx_order,dx,X,
-                                   len,deriv_var,dep_var_idx,indep_var_idx)
+                                   len_of_indep_vars,deriv_var,dep_var_idx,indep_var_idx)
                 return :(broadcast($name, $aux))
             else
                 aux_1 = discretize_2(input.args[1],deriv_order,approx_order,dx,X,
-                                     len,deriv_var,dep_var_idx,indep_var_idx)
+                                     len_of_indep_vars,deriv_var,dep_var_idx,indep_var_idx)
                 aux_2 = discretize_2(input.args[2],deriv_order,approx_order,dx,X,
-                                     len,deriv_var,dep_var_idx,indep_var_idx)
+                                     len_of_indep_vars,deriv_var,dep_var_idx,indep_var_idx)
                 return :(broadcast($name, $aux_1, $aux_2))
             end
         end
@@ -130,7 +230,7 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
     domain = Array{Any}(undef,num_indep_vars)
     dx = Array{Any}(undef,num_indep_vars)
     X = Array{Any}(undef,num_indep_vars)
-    len = Array{Any}(undef,num_indep_vars)
+    len_of_indep_vars = Array{Any}(undef,num_indep_vars)
     k = 0
     for i = 1:num_indep_vars
         var = nameof(pdesys.domain[i].variables)
@@ -139,11 +239,11 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
         if var != :(t)
             dx[i] = discretization.dxs[i-k]
             X[i] = domain[i].lower:dx[i]:domain[i].upper
-            len[i] = size(X[i],1)
+            len_of_indep_vars[i] = size(X[i],1)
         else
             dx[i] = 0.0
             X[i] = 0.0
-            len[i] = 0.0
+            len_of_indep_vars[i] = 0.0
             tdomain = pdesys.domain[1].domain
             @assert tdomain isa IntervalDomain
             k = 1
@@ -178,48 +278,51 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
         dep_var_idx[var] = j
     end
     for (var,j) in dep_var_idx
-        aux = discretize_2( eqs[j].rhs,0,approx_order,dx,X,len,
+        aux = discretize_2( eqs[j].rhs,0,approx_order,dx,X,len_of_indep_vars,
                             [],dep_var_idx,indep_var_idx)
-        dep_var_disc[var] = @RuntimeGeneratedFunction(:((all_bcs,u,t) -> $aux))
+        dep_var_disc[var] = @RuntimeGeneratedFunction(:((Q,u,t) -> $aux))
     end
 
     ### Declare and define boundary conditions ################################
 
     # TODO: extend to Neumann BCs and Robin BCs
     lhs_deriv_depvars_bcs = get_bcs(pdesys.bcs,tdomain,domain[2])
-    t = 0.0
-    u_ic = Array{Float64}(undef,len[2]-2,num_dep_vars)
-    u_left_bc = Array{Any}(undef,num_dep_vars)
-    u_right_bc = Array{Any}(undef,num_dep_vars)
-    all_bcs = Array{RobinBC}(undef,num_dep_vars)
+    u_ic = Array{Float64}(undef,len_of_indep_vars[2]-2,num_dep_vars)
+    robin_bc_func = Array{Any}(undef,num_dep_vars)
+    Q = Array{RobinBC}(undef,num_dep_vars)
 
     for var in keys(dep_var_idx)
         j = dep_var_idx[var]
         bcs = lhs_deriv_depvars_bcs[var]
 
-        ic = @RuntimeGeneratedFunction(:((x,t) -> $(bcs[1])))
-        u_ic[:,j] = ic.(X[2][2:len[2]-1],t)
+        # Initial condition depends on space but not time
+        ic = @RuntimeGeneratedFunction(:(x -> $(bcs["ic"])))
+        u_ic[:,j] = ic.(X[2][2:len_of_indep_vars[2]-1])
 
-        left_bc_fn = :((x, t) -> $(bcs[2]))
-        right_bc_fn = :((x, t) -> $(bcs[3]))
-        u_left_bc[j] = @RuntimeGeneratedFunction(left_bc_fn)
-        u_right_bc[j] = @RuntimeGeneratedFunction(right_bc_fn)
+        # Boundary conditions depend on time and so will be evaluated at t within the
+        # ODE function for the discretized PDE (below)
+        # We use @RuntimeGeneratedFunction to do this efficiently
+        αl, βl, γl = bcs["left bc"]
+        αr, βr, γr = bcs["right bc"]
+        # Right now there is only one independent variable, so dx is always dx[2]
+        # i.e. the spatial variable
+        robin_bc_func[j] = @RuntimeGeneratedFunction(:(t -> begin
+            RobinBC(($(αl), $(βl), $(γl)), ($(αr), $(βr), $(γr)), $(dx[2]))
+        end))
     end
 
     ### Define the discretized PDE as an ODE function #########################
 
     function f(du,u,p,t)
 
-        # Boundary conditions can vary with respect to time
+        # Boundary conditions can vary with respect to time (but not space)
         for j in 1:num_dep_vars
-            left_bc_eval = u_left_bc[j](X[2][1],0.0)
-            right_bc_eval = u_right_bc[j](last(X[2]),0.0)
-            all_bcs[j] = DirichletBC(left_bc_eval,right_bc_eval)
+            Q[j] = robin_bc_func[j](t)
         end
 
         for (var,disc) in dep_var_disc
             j = dep_var_idx[var]
-            res = disc(all_bcs,u,t)
+            res = disc(Q,u,t)
             if haskey(lhs_deriv_depvars,var)
                 du[:,j] = res
             else
@@ -230,5 +333,13 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
     end
 
     # Return problem ##########################################################
-    return PDEProblem(ODEProblem(f,u_ic,(tdomain.lower,tdomain.upper),nothing),all_bcs,X)
+    # The second entry, robin_bc_func, is stored as the "extrapolation" in the
+    # PDEProblem. Since this can in general be time-dependent, it must be evaluated
+    # at the correct time when used, e.g.
+    #   prob.extrapolation[1](t[i])*sol[:,1,i]
+    return PDEProblem(
+        ODEProblem(f,u_ic,(tdomain.lower,tdomain.upper),nothing),
+        robin_bc_func,
+        X
+    )
 end
diff --git a/src/exceptions.jl b/src/exceptions.jl
new file mode 100644
index 000000000..06d5fc82c
--- /dev/null
+++ b/src/exceptions.jl
@@ -0,0 +1,7 @@
+struct BoundaryConditionError <: Exception
+    msg::String
+end
+
+Base.showerror(io::IO, e::BoundaryConditionError) = print(
+    io, "BoundaryConditionError: ", e.msg
+)
\ No newline at end of file
diff --git a/test/MOL_1D_Linear_Diffusion.jl b/test/MOL_1D_Linear_Diffusion.jl
index 9f2f300fc..22a72f77f 100644
--- a/test/MOL_1D_Linear_Diffusion.jl
+++ b/test/MOL_1D_Linear_Diffusion.jl
@@ -3,10 +3,13 @@
 # TODO: Add more complex tests.
 
 # Packages and inclusions
-using ModelingToolkit,DiffEqOperators,DiffEqBase,LinearAlgebra,Test
+using ModelingToolkit,DiffEqOperators,DiffEqBase,LinearAlgebra,Test,OrdinaryDiffEq
 
 # Tests
 @testset "Test 00: Dt(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(..)
@@ -15,13 +18,13 @@ using ModelingToolkit,DiffEqOperators,DiffEqBase,LinearAlgebra,Test
 
     # 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,Float64(pi)) ~ -exp(-t)]
 
     # Space and time domains
     domains = [t ∈ IntervalDomain(0.0,1.0),
-               x ∈ IntervalDomain(0.0,1.0)]
+               x ∈ IntervalDomain(0.0,Float64(pi))]
 
     # PDE system
     pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
@@ -35,22 +38,24 @@ using ModelingToolkit,DiffEqOperators,DiffEqBase,LinearAlgebra,Test
     prob = discretize(pdesys,discretization)
 
     # Solve ODE problem
-    using OrdinaryDiffEq
     sol = solve(prob,Tsit5(),saveat=0.1)
+    x = prob.space[2]
+    t = sol.t
 
     # Plot and save results
     # using Plots
-    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
+    # plot()
+    # for i in 1:4
+    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
+    #     scatter!(x, u_exact(x, t[i]))
+    # end
     # savefig("MOL_1D_Linear_Diffusion_Test00.png")
 
-    # Test
-    n = size(sol,1)
-    t_f = size(sol,3)
-
-    @test sol[:,1,t_f] ≈ zeros(n) atol = 0.001;
+    # Test against exact solution
+    for i in 1:size(t,1)
+        u_approx = Array(prob.extrapolation[1](t[i])*sol.u[i])
+        @test u_approx ≈ u_exact(x, t[i]) atol=0.01
+    end
 end
 
 @testset "Test 01: Dt(u(t,x)) ~ D*Dxx(u(t,x))" begin
@@ -84,15 +89,17 @@ end
     prob = discretize(pdesys,discretization)
 
     # Solve ODE problem
-    using OrdinaryDiffEq
     sol = solve(prob,Tsit5(),saveat=0.1)
 
     # Plot and save results
+    # x = prob.space[2]
+    # t = sol.t
+    
     # using Plots
-    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
+    # plot()
+    # for i in 1:4
+    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
+    # end
     # savefig("MOL_1D_Linear_Diffusion_Test01.png")
 
     # Test
@@ -137,19 +144,331 @@ end
     prob = discretize(pdesys,discretization)
 
     # Solve ODE problem
-    using OrdinaryDiffEq
     sol = solve(prob,Tsit5(),saveat=0.1)
     
     # Plot and save results
+    # x = prob.space[2]
+    # t = sol.t
+    
     # using Plots
-    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
+    # plot()
+    # for i in 1:4
+    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
+    # end
     # savefig("MOL_1D_Linear_Diffusion_Test02.png")
 
     # Test
     n = size(sol,1)
     t_f = size(sol,3)
-    @test sol[:,1,t_f] ≈ zeros(n) atol = 0.1;
+    @test sol[:,1,t_f] ≈ zeros(n) atol=0.01;
+end
+
+@testset "Test 03: Dt(u(t,x)) ~ Dxx(u(t,x)), Neumann BCs" begin
+    # 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),
+           Dx(u(t,0)) ~ exp(-t),
+           Dx(u(t,Float64(pi))) ~ -exp(-t)]
+
+    # Space and time domains
+    domains = [t ∈ IntervalDomain(0.0,1.0),
+               x ∈ IntervalDomain(0.0,Float64(pi))]
+
+    # PDE system
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+
+    # Method of lines discretization
+    dx = 0.1
+    order = 2
+    discretization = MOLFiniteDifference(dx,order)
+
+    # Convert the PDE problem into an ODE problem
+    prob = discretize(pdesys,discretization)
+
+    # Solve ODE problem
+    sol = solve(prob,Tsit5(),saveat=0.1)
+    x = prob.space[2]
+    t = sol.t
+
+    # Plot and save results
+    # using Plots
+    # plot()
+    # for i in 1:4
+    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
+    #     scatter!(x, u_exact(x, t[i]))
+    # end
+    # savefig("MOL_1D_Linear_Diffusion_Test03.png")
+
+    # Test against exact solution
+    for i in 1:size(t,1)
+        u_approx = Array(prob.extrapolation[1](t[i])*sol.u[i])
+        @test u_approx ≈ u_exact(x, t[i]) atol=0.01
+    end
+end
+
+@testset "Test 04: Dt(u(t,x)) ~ Dxx(u(t,x)), Neumann + Dirichlet BCs" begin
+    # 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,0) ~ 0.0,
+           Dx(u(t,Float64(pi))) ~ -exp(-t)]
+
+    # Space and time domains
+    domains = [t ∈ IntervalDomain(0.0,1.0),
+               x ∈ IntervalDomain(0.0,Float64(pi))]
+
+    # PDE system
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+
+    # Method of lines discretization
+    dx = 0.1
+    order = 2
+    discretization = MOLFiniteDifference(dx,order)
+
+    # Convert the PDE problem into an ODE problem
+    prob = discretize(pdesys,discretization)
+
+    # Solve ODE problem
+    sol = solve(prob,Tsit5(),saveat=0.1)
+    x = prob.space[2]
+    t = sol.t
+
+    # Plot and save results
+    # using Plots
+    # plot()
+    # for i in 1:4
+    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
+    #     scatter!(x, u_exact(x, t[i]))
+    # end
+    # savefig("MOL_1D_Linear_Diffusion_Test04.png")
+
+    # Test against exact solution
+    for i in 1:size(t,1)
+        u_approx = Array(prob.extrapolation[1](t[i])*sol.u[i])
+        @test u_approx ≈ u_exact(x, t[i]) atol=0.01
+    end
+end
+
+@testset "Test 05: Dt(u(t,x)) ~ Dxx(u(t,x)), Robin BCs" begin
+    # 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)),
+           4u(t,1.0) + Dx(u(t,1.0)) ~ exp(-t) * (4sin(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
+    dx = 0.01
+    order = 2
+    discretization = MOLFiniteDifference(dx,order)
+
+    # Convert the PDE problem into an ODE problem
+    prob = discretize(pdesys,discretization)
+
+    # Solve ODE problem
+    sol = solve(prob,Tsit5(),saveat=0.1)
+    x = prob.space[2]
+    t = sol.t
+
+    # Plot and save results
+    # using Plots
+    # plot()
+    # for i in 1:4
+    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
+    #     scatter!(x, u_exact(x, t[i]))
+    # end
+    # savefig("MOL_1D_Linear_Diffusion_Test05.png")
+
+    # Test against exact solution
+    for i in 1:size(t,1)
+        u_approx = Array(prob.extrapolation[1](t[i])*sol.u[i])
+        @test u_approx ≈ u_exact(x, t[i]) atol=0.05
+    end
+end
+
+@testset "Test 06: Dt(u(t,x)) ~ Dxx(u(t,x)), time-dependent Robin BCs" begin
+    # 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),
+           t^2 * u(t,-1.0) + 3Dx(u(t,-1.0)) ~ exp(-t) * (t^2 * sin(-1.0) + 3cos(-1.0)),
+           4u(t,1.0) + t * Dx(u(t,1.0)) ~ exp(-t) * (4sin(1.0) + t * 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
+    dx = 0.01
+    order = 2
+    discretization = MOLFiniteDifference(dx,order)
+
+    # Convert the PDE problem into an ODE problem
+    prob = discretize(pdesys,discretization)
+
+    # Solve ODE problem
+    sol = solve(prob,Tsit5(),saveat=0.1)
+    x = prob.space[2]
+    t = sol.t
+
+    # Plot and save results
+    # using Plots
+    # plot()
+    # for i in 1:4
+    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
+    #     scatter!(x, u_exact(x, t[i]))
+    # end
+    # savefig("MOL_1D_Linear_Diffusion_Test06.png")
+
+    # Test against exact solution
+    for i in 1:size(t,1)
+        u_approx = Array(prob.extrapolation[1](t[i])*sol.u[i])
+        @test u_approx ≈ u_exact(x, t[i]) atol=0.05
+    end
+end
+
+@testset "Test errors" begin
+    # Parameters, variables, and derivatives
+    @parameters t x
+    @variables u(..) v(..)
+    @derivatives Dt'~t
+    @derivatives Dx'~x
+    @derivatives Dxx''~x
+
+    # 1D PDE and boundary conditions
+    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
+    # Space and time domains
+    domains = [t ∈ IntervalDomain(0.0,1.0), x ∈ IntervalDomain(0.0,1.0)]
+
+    # Method of lines discretization
+    dx = 0.1
+    order = 2
+    discretization = MOLFiniteDifference(dx,order)
+
+    # Missing boundary condition
+    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
+           Dx(u(t,1)) ~ 0.0]    
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+    @test_throws BoundaryConditionError discretize(pdesys,discretization)
+    
+    # Boundary condition not at t=0
+    bcs = [u(1,x) ~ 0.5 + sin(2pi*x),
+           Dx(u(t,0)) ~ 0.0,
+           Dx(u(t,1)) ~ 0.0]
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+    @test_throws BoundaryConditionError discretize(pdesys,discretization)
+
+    # Boundary condition not at an edge of the domain
+    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
+           Dx(u(t,0)) ~ 0.0,
+           Dx(u(t,0.5)) ~ 0.0]    
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+    @test_throws BoundaryConditionError discretize(pdesys,discretization)
+    
+    # Second-order derivative in BC
+    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
+           Dxx(u(t,0)) ~ 0.0,
+           Dx(u(t,0.5)) ~ 0.0]    
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+    @test_throws BoundaryConditionError discretize(pdesys,discretization)
+   
+    # Wrong format for Robin BCs
+    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
+           Dx(u(t,0)) ~ 0.0,
+           u(t,1) * Dx(u(t,1)) ~ 0.0]    
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+    @test_throws BoundaryConditionError discretize(pdesys,discretization)
+    
+    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
+           Dx(u(t,0)) ~ 0.0,
+           u(t,1) + Dxx(u(t,1)) ~ 0.0]    
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+    @test_throws BoundaryConditionError discretize(pdesys,discretization)
+    
+    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
+           Dx(u(t,0)) ~ 0.0,
+           u(t,1) + 2Dxx(u(t,1)) ~ 0.0]    
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+    @test_throws BoundaryConditionError discretize(pdesys,discretization)
+    
+    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
+           Dx(u(t,0)) ~ 0.0,
+           Dx(u(t,1)) + u(t,1) ~ 0.0]    
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+    @test_throws BoundaryConditionError discretize(pdesys,discretization)
+    
+    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
+           Dx(u(t,0)) ~ 0.0,
+           u(t,1) / 2 + Dx(u(t,1)) ~ 0.0]    
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+    @test_throws BoundaryConditionError discretize(pdesys,discretization)
+    
+    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
+           Dx(u(t,0)) ~ 0.0,
+           u(t,1) + Dx(u(t,1)) / 2 ~ 0.0]    
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+    @test_throws BoundaryConditionError discretize(pdesys,discretization)
+    
+    # Mismatching arguments
+    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
+    Dx(u(t,0)) ~ 0.0,
+    u(t,0) + Dx(u(t,1)) ~ 0.0]    
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+    @test_throws BoundaryConditionError discretize(pdesys,discretization)
+    
+    # Mismatching variables
+    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
+    Dx(u(t,0)) ~ 0.0,
+    u(t,1) + Dx(v(t,1)) ~ 0.0]    
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+    @test_throws BoundaryConditionError discretize(pdesys,discretization)
+    
 end