diff --git a/REQUIRE b/REQUIRE
index 02af2de..d14e867 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -2,3 +2,4 @@ julia 0.3
 Compat
 Docile
 ODE 0.2.1
+FastAnonymous
diff --git a/src/QuDynamics.jl b/src/QuDynamics.jl
index a668def..5da7617 100644
--- a/src/QuDynamics.jl
+++ b/src/QuDynamics.jl
@@ -4,6 +4,7 @@ module QuDynamics
     using ODE
     using Expokit
     using ExpmV
+    using FastAnonymous
     VERSION < v"0.4-" && using Docile
     include("quequations.jl")
     include("propmachinery.jl")
diff --git a/src/propodesolvers.jl b/src/propodesolvers.jl
index 4664ccb..03719bb 100644
--- a/src/propodesolvers.jl
+++ b/src/propodesolvers.jl
@@ -31,7 +31,8 @@ for (qu_ode_type,ode_solver) in type_to_method_ode
             dims = size(current_qustate)
             # Convert the current_qustate to complex as it might result in a Inexact Error. After complex is in QuBase.jl (PR #38)
             # we could just do a complex(vec(current_qustate)) avoiding the coeffs(coeffs(vec(current_qustate))).
-            next_state = $ode_solver((t,y)-> -im*coeffs(op)*y, complex(coeffs(vec(current_qustate))), [current_t, t], points=:specified,
+            f = @anon (t,y)-> scale!(coeffs(op)*y, -im)
+            next_state = $ode_solver(f, complex(coeffs(vec(current_qustate))), [current_t, t], points=:specified,
                                   reltol = get(prob.options, :reltol, 1.0e-5), abstol = get(prob.options, :abstol, 1.0e-8))[2][end]
             CQST = QuBase.similar_type(current_qustate)
             return CQST(reshape(next_state, dims), bases(current_qustate))