diff --git a/src/MOL_discretization.jl b/src/MOL_discretization.jl
index 25c508697..50cd63e40 100644
--- a/src/MOL_discretization.jl
+++ b/src/MOL_discretization.jl
@@ -1,4 +1,3 @@
-using Printf
 using ModelingToolkit: operation
 # Method of lines discretization scheme
 struct MOLFiniteDifference{T} <: DiffEqBase.AbstractDiscretization
@@ -7,6 +6,11 @@ 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 not recognize 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) = γ
@@ -16,7 +20,6 @@ end
 # Dx(u(t, 0)) ~ γ                       ---> (0, 1, γ)
 # α * u(t, 0)) + β * Dx(u(t, 0)) ~ γ    ---> (α, β, γ)
 # In general, all of α, β, and γ could be Expr 
-# For now, we restrict to the case where α and β are Float64
 function get_bcs(bcs,tdomain,domain)
     lhs_deriv_depvars_bcs = Dict()
     num_bcs = size(bcs,1)
@@ -35,8 +38,47 @@ function get_bcs(bcs,tdomain,domain)
             α = 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 = operation(lhs_l).name
+                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
+                var_l = operation(lhs_l.args[2]).name
+                bc_args_l = lhs_l.args[2].args
+            else
+                throw_bc_err(bcs[i])
+            end
+            # Right side of the expression should be Differential or β * Differential
+            if operation(lhs_r) isa Differential
+                β = 1.0
+                var_r = operation(lhs_r.args[1]).name
+                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(β) : β
+                @assert operation(lhs_r.args[2]) isa Differential
+                var_r = operation(lhs_r.args[2].args[1]).name
+                bc_args_r = lhs_r.args[2].args[1].args
+            else
+                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
+            var = var_l
+            @assert bc_args_l == bc_args_r
+            bc_args = bc_args_l
         else
-            # TODO: RobinBC
+            throw_bc_err(bcs[i])
         end
         if !haskey(lhs_deriv_depvars_bcs,var)
             # Initialize dict of boundary conditions for this variable
@@ -52,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=0, x=$(domain.lower), or x=$(domain.upper)"
             ))
         end
         # Create value
@@ -70,7 +112,7 @@ function get_bcs(bcs,tdomain,domain)
     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(@sprintf "missing %s for %s" key var))
+                throw(BoundaryConditionError("missing $key for $var"))
             end
         end
     end
diff --git a/test/MOL_1D_Linear_Diffusion.jl b/test/MOL_1D_Linear_Diffusion.jl
index 4ff6047be..948a4816a 100644
--- a/test/MOL_1D_Linear_Diffusion.jl
+++ b/test/MOL_1D_Linear_Diffusion.jl
@@ -163,7 +163,7 @@ end
     @test sol[:,1,t_f] ≈ zeros(n) atol = 0.1;
 end
 
-@testset "Test 03: Dt(u(t,x)) ~ Dxx(u(t,x)), Dx(u(t,0)) ~ 0, Dx(u(t,1)) ~ 0" begin
+@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)
 
@@ -200,6 +200,112 @@ end
     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)),
+           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
+    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()
@@ -207,19 +313,19 @@ end
         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")
+    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.01
+        @test u_approx ≈ u_exact(x, t[i]) atol=0.1
     end
 end
 
 @testset "Test errors" begin
     # Parameters, variables, and derivatives
     @parameters t x
-    @variables u(..)
+    @variables u(..) v(..)
     @derivatives Dt'~t
     @derivatives Dx'~x
     @derivatives Dxx''~x
@@ -234,6 +340,12 @@ end
     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,
@@ -247,10 +359,44 @@ end
            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)
     
-    # Missing boundary condition
     bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
-           Dx(u(t,1)) ~ 0.0]    
+           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 AssertionError 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 AssertionError discretize(pdesys,discretization)
+    
 end