diff --git a/Project.toml b/Project.toml
index 0e9ff2f1..1d174520 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,6 +14,7 @@ IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 JuAFEM = "30d91d44-8115-11e8-1d28-c19a5ac16de8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
@@ -22,7 +23,6 @@ SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8"
 Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
 VoronoiDelaunay = "72f80fcb-8c52-57d9-aff0-40c1a3526986"
 
@@ -35,12 +35,12 @@ GeometricalPredicates = "0.3, 0.4"
 Interpolations = "0.11, 0.12"
 IterativeSolvers = "0.7, 0.8"
 LinearMaps = "2"
+ModelingToolkit  = "3.0.1"
 NearestNeighbors = "0.4.4"
 OrdinaryDiffEq = "5.18"
 ProgressMeter = "1"
 RecipesBase = "0.7, 0.8, 1.0"
 StaticArrays = "0.10, 0.11, 0.12"
-SymEngine = "0.6, 0.7, 0.8"
 Tensors = "1.0.1"
 VoronoiDelaunay = "0.4"
 julia = "1"
diff --git a/src/CoherentStructures.jl b/src/CoherentStructures.jl
index ef0fae44..97e153a2 100644
--- a/src/CoherentStructures.jl
+++ b/src/CoherentStructures.jl
@@ -41,7 +41,9 @@ const VD = VoronoiDelaunay
 import JuAFEM
 const JFM = JuAFEM
 using RecipesBase
-import SymEngine
+import ModelingToolkit
+using ModelingToolkit: Variable, Differential, simplified_expr,
+      expand_derivatives, Expression, Operation, Constant
 
 # contains a list of exported functions
 include("exports.jl")
diff --git a/src/streammacros.jl b/src/streammacros.jl
index 84551493..cb154335 100644
--- a/src/streammacros.jl
+++ b/src/streammacros.jl
@@ -196,7 +196,7 @@ end
 function streamline_derivatives(H::Symbol, formulas::Expr)
     # symbols that are not supposed to be substituted
     # (additional to symbols defined in Base)
-    bound_symbols = [:x, :y, :p, :t, keys(diff_dict)...]
+    bound_symbols = [:x, :y, :p, :t]
     H = substitutions(H, formulas, bound_symbols)
 
     # symbolic gradient and hessian (note the broadcast)
@@ -216,90 +216,32 @@ function streamline_derivatives(H::Symbol, formulas::Expr)
 end
 
 ########################################################################################
-#            symbolic differentiation of expressions using SymEngine                   #
+#            symbolic differentiation of expressions using ModelingToolkit             #
 ########################################################################################
 
-sgn(x) = (x > 0) ? one(x) : (x < 0) ? -one(x) : zero(x)
-heaviside(x) = 0 < x ? one(x) : zero(x)
-window(x, a, b) = (a <= x < b) ? one(x) : zero(x)
-
-# manually define  derivatives for functions that SymEngine cant differentiate
-diff_dict = Dict()
-diff_dict[:abs] = :sgn
-diff_dict[:sgn] = :zero
-diff_dict[:heaviside] = :zero
-diff_dict[:zero] = :zero
-diff_dict[:window] = :zero
+ModelingToolkit.derivative(::typeof(sign), x, ::Val{1}) = 0
 
 function expr_diff(expr::Expr, var::Symbol)
-    # not a nice way to differentiate expressions, but ReverseDiffSource
-    # is broken.
-    expr_sym = SymEngine.Basic(expr)
-    d_expr_sym = SymEngine.diff(expr_sym, var)
-    d_expr = Meta.parse(SymEngine.toString.(d_expr_sym))
-
-    # resolve derivatives that SymEngine doesn't know using diff_dict
-    d_expr = additional_derivatives(d_expr)
+    D = Differential(var) 
 
-    # clean up unresolved substitutions that result from SymEnginge treating
-    # unknown derivatives
-    d_expr = substitution_cleaner(d_expr)
+    toolkit_expr = convert(Expression, expr) 
+    d_expr = simplified_expr(expand_derivatives(D(toolkit_expr)))
+    d_expr = convert(Expression, simple_simplifier(d_expr))
+    return simplified_expr(propagate_constants(d_expr))
 
-    # clean up zeros
-    d_expr = simple_simplifier(d_expr)
 end
 expr_diff(expr, var::Symbol) = expr == var ? 1 : 0
 
-function additional_derivatives(expr::Expr)
-    # some functions like abs(x) are not treated by SymEngine and block the
-    # expression. For example, diff(Basic(:(abs(x^2+1)),:x)) returns a SymEngine Object
-    # whose string representation is parsed to :(Derivative(abs(1 + x ^ 2), x)),
-    # whose AST is:
-    # Expr
-    #  head: Symbol call
-    #  args: Array{Any}((3,))
-    #    1: Symbol Derivative
-    #    2: Expr
-    #           ... Body of expression ...
-    #    3: Symbol x
-    # typ: Any
-
-    # detect expressions of this form
-    if expr.head === :call && expr.args[1] === :Derivative
-        f = expr.args[2].args[1]
-        var = expr.args[3]
-        f_arg = expr.args[2].args[2]
-        if haskey(diff_dict, f)  # try if diff_dict provides a rescue
-            df = diff_dict[f]
-            inner_d = expr_diff(f_arg, var)
-            df_computed_manually = :($df($f_arg) * $inner_d)
-            return additional_derivatives(df_computed_manually)       # handle nested problems
-        end
-    end
-    return Expr(expr.head, additional_derivatives.(expr.args)...) # call recursively on subexpressions
-end
-additional_derivatives(expr) = expr
-
-# A second thing that SymEngine does is returning expressions of the form
-# Subs(ex1, symb1, ex2). Resolve these substitutions
-function substitution_cleaner(expr::Expr)
-    if expr.head === :call && expr.args[1] === :Subs
-        return substitution_cleaner(sym_subst(
-            expr.args[2],
-            expr.args[3],
-            expr.args[4],
-        ))
-    end
-    return Expr(expr.head, substitution_cleaner.(expr.args)...)
-end
-substitution_cleaner(expr) = expr
-
 # perform some basic simplifications like getting rid of ones and zeros
 function simple_simplifier(expr::Expr)
     args = simple_simplifier.(expr.args)
     if expr.head !== :call
         return Expr(expr.head, args...)
     end
+    if args[1] === :one
+        return 1
+    end
+
     if args[1] === :zero
         return 0
     end
@@ -333,10 +275,22 @@ function simple_simplifier(expr::Expr)
     end
     return Expr(expr.head, args...)
 end
-
 simple_simplifier(expr) = expr
-expr_grad(expr, coord_vars::Vector{Symbol}) = expr_diff.(expr, coord_vars)
 
+function propagate_constants(expr::ModelingToolkit.Operation)
+        if all(isconstant.(expr.args))
+            return expr.op(propagate_constants.(expr.args)...) 
+        else
+            return Operation(expr.op, propagate_constants.(expr.args))
+        end
+end
+propagate_constants(x::Constant) = x.value
+propagate_constants(x) = x
+
+isconstant(x::Operation) = !(x.op isa Variable) && all(isconstant.(x.args))
+isconstant(x::Constant) = true
+
+expr_grad(expr, coord_vars::Vector{Symbol}) = expr_diff.(expr, coord_vars)
 function hessian(expr, coord_vars)
     ∇expr = expr_grad(expr, coord_vars)
     ∇²expr = expr_grad(expr)
diff --git a/src/velocityfields.jl b/src/velocityfields.jl
index db96e5d1..e05bafbd 100644
--- a/src/velocityfields.jl
+++ b/src/velocityfields.jl
@@ -26,6 +26,7 @@ bickleyJetEqVari! = ODE.ODEFunction{true}((DU, U, p, t) -> DU .= bickleyJetEqVar
 # rotating double gyre flow  [Mosovsky & Meiss, 2011]
 @define_stream Ψ_rot_dgyre begin
     st          = heaviside(t)*heaviside(1-t)*t^2*(3-2*t) + heaviside(t-1)
+    heaviside(x)= 0.5*(sign(x) + 1)
     Ψ_P         = sin(2π*x)*sin(π*y)
     Ψ_F         = sin(π*x)*sin(2π*y)
     Ψ_rot_dgyre = (1-st) * Ψ_P + st * Ψ_F