diff --git a/Project.toml b/Project.toml
index e1d96a9f5..18bd4434c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "DataDrivenDiffEq"
 uuid = "2445eb08-9709-466a-b3fc-47e12bd697a2"
 authors = ["Julius Martensen <julius.martensen@gmail.com>"]
-version = "0.6.0"
+version = "0.6.1"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
diff --git a/docs/src/prob_and_solve.md b/docs/src/prob_and_solve.md
index 73b23720f..915ba4fd0 100644
--- a/docs/src/prob_and_solve.md
+++ b/docs/src/prob_and_solve.md
@@ -67,6 +67,15 @@ ps = parameter_map(res)
 
 ## [Optional Arguments](@id optional_arguments)
 
+!!! info
+The keyword argument `eval_expression` controls the function creation
+behavior. `eval_expression=true` means that `eval` is used, so normal
+world-age behavior applies (i.e. the functions cannot be called from
+the function that generates them). If `eval_expression=false`,
+then construction via GeneralizedGenerated.jl is utilized to allow for
+same world-age evaluation. However, this can cause Julia to segfault
+on sufficiently large basis functions. By default eval_expression=false.
+
 Koopman based algorithms can be called without a [`Basis`](@ref), resulting in dynamic mode decomposition like methods, or with a basis for extened dynamic mode decomposition :
 
 ```julia
@@ -81,6 +90,9 @@ Possible keyworded arguments include
 + `digits` controls the digits / rounding used for deriving the system equations (`digits = 1` would round `10.02` to `10.0`)
 + `operator_only` returns a `NamedTuple` containing the operator, input and output mapping and matrices used for updating the operator as described [here](https://arxiv.org/pdf/1406.7187.pdf)
 
+!!! info
+If `eval_expression` is set to `true`, the returning result of the Koopman based inference will not contain a parametrized equation, but rather use the numeric values of the operator/generator.
+
 SINDy based algorithms can be called like :
 
 ```julia
diff --git a/src/optimizers/sparseregression.jl b/src/optimizers/sparseregression.jl
index e86266708..fc7bca917 100644
--- a/src/optimizers/sparseregression.jl
+++ b/src/optimizers/sparseregression.jl
@@ -3,11 +3,17 @@
 $(SIGNATURES)
 
 Implements a sparse regression, given an `AbstractOptimizer` or `AbstractSubspaceOptimizer`.
-`maxiter` indicate the maximum iterations for each call of the optimizer, `abstol` the absolute tolerance of
+`X` denotes the coefficient matrix, `A` the design matrix and `Y` the matrix of observed or target values.
+`X` can be derived via `init(opt, A, Y)`.
+`maxiter` indicates the maximum iterations for each call of the optimizer, `abstol` the absolute tolerance of
 the difference between iterations in the 2 norm. If the optimizer is called with a `Vector` of thresholds, each `maxiter` indicates
 the maximum iterations for each threshold.
 
-If `progress` is set to `true`, a progressbar will be available.
+If `progress` is set to `true`, a progressbar will be available. `progress_outer` and `progress_offset` are used to compute the initial offset of the
+progressbar.
+
+If used with a `Vector` of thresholds, the functions `f` with signature `f(X, A, Y)` and `g` with signature `g(x, threshold) = G(f(X, A, Y))` with the arguments given as stated above can be passed in. These are
+used for finding the pareto-optimal solution to the sparse regression. 
 """
 function sparse_regression!(X, A, Y, opt::AbstractOptimizer{T};
     maxiter::Int = maximum(size(A)),
@@ -83,7 +89,7 @@ function sparse_regression!(X, A, Y, opt::AbstractOptimizer{T};
             )
         end
     end
-    
+
     for i in 1:size(Y, 2)
         @views clip_by_threshold!(X[:, i], λs[i])
     end
diff --git a/src/solution.jl b/src/solution.jl
index 7d927c832..ea696b1d1 100644
--- a/src/solution.jl
+++ b/src/solution.jl
@@ -108,7 +108,7 @@ function Base.print(io::IO, r::DataDrivenSolution, fullview::DataType)
     println(io, "")
     print(io, r.res)
     println(io, "")
-    if length(r.res.ps) > 0 
+    if length(r.res.ps) > 0
         x = parameter_map(r)
         println(io, "Parameters:")
         for v in x
@@ -154,14 +154,15 @@ end
 
 
 # Explicit sindy
-function build_solution(prob::DataDrivenProblem, Ξ::AbstractMatrix, opt::Optimize.AbstractOptimizer, b::Basis)
+function build_solution(prob::DataDrivenProblem, Ξ::AbstractMatrix, opt::Optimize.AbstractOptimizer, b::Basis;
+    eval_expression = false)
     if all(iszero(Ξ))
         @warn "Sparse regression failed! All coefficients are zero."
         return DataDrivenSolution(
-        nothing , :failed, nothing, opt, Ξ, (Problem = prob, Basis = b, nothing), 
+        nothing , :failed, nothing, opt, Ξ, (Problem = prob, Basis = b, nothing),
     )
     end
-    
+
     eqs, ps, p_ = build_parametrized_eqs(Ξ, b)
 
     # Build the lhs
@@ -176,7 +177,8 @@ function build_solution(prob::DataDrivenProblem, Ξ::AbstractMatrix, opt::Optimi
         eqs, states(b),
         parameters = [parameters(b); p_], iv = independent_variable(b),
         controls = controls(b), observed = observed(b),
-        name = gensym(:Basis)
+        name = gensym(:Basis),
+        eval_expression = eval_expression
     )
 
     sparsity = norm(Ξ, 0)
@@ -225,15 +227,15 @@ function build_solution(prob::DataDrivenProblem, Ξ::AbstractMatrix, opt::Optimi
 end
 
 function build_solution(prob::DataDrivenProblem, Ξ::AbstractMatrix, opt::Optimize.AbstractSubspaceOptimizer,
-    b::Basis, implicits::Vector{Num})
+    b::Basis, implicits::Vector{Num}; eval_expression = false)
 
     if all(iszero(Ξ))
         @warn "Sparse regression failed! All coefficients are zero."
         return DataDrivenSolution(
-        nothing , :failed, nothing, opt, Ξ, (Problem = prob, Basis = b, nothing), 
+        nothing , :failed, nothing, opt, Ξ, (Problem = prob, Basis = b, nothing),
     )
     end
-    
+
     eqs, ps, p_ = build_parametrized_eqs(Ξ, b)
     eqs = [0 .~ eq for eq in eqs]
 
@@ -242,7 +244,8 @@ function build_solution(prob::DataDrivenProblem, Ξ::AbstractMatrix, opt::Optimi
         eqs, states(b),
         parameters = [parameters(b); p_], iv = independent_variable(b),
         controls = controls(b), observed = observed(b),
-        name = gensym(:Basis)
+        name = gensym(:Basis),
+        eval_expression = eval_expression
     )
 
     sparsity = norm(Ξ, 0)
@@ -298,8 +301,8 @@ function _round!(x::AbstractArray{T, N}, digits::Int) where {T, N}
 end
 
 #function build_solution(prob::DataDrivenProblem, Ξ::AbstractMatrix, opt::AbstractKoopmanAlgorithm)
-function build_solution(prob::DataDrivenProblem, k, C, B, Q, P, inds, b::AbstractBasis, alg::AbstractKoopmanAlgorithm; digits::Int = 10)
-    # Build parameterized equations
+function build_solution(prob::DataDrivenProblem, k, C, B, Q, P, inds, b::AbstractBasis, alg::AbstractKoopmanAlgorithm; digits::Int = 10, eval_expression = false)
+    # Build parameterized equations, inds indicate the location of basis elements containing an input
     Ξ = zeros(eltype(B), size(C,2), length(b))
 
     Ξ[:, inds] .= real.(Matrix(k))
@@ -308,7 +311,14 @@ function build_solution(prob::DataDrivenProblem, k, C, B, Q, P, inds, b::Abstrac
     end
 
     # Transpose because of the nature of build_parametrized_eqs
-    eqs, ps, p_ = build_parametrized_eqs(_round!(C*Ξ, digits)', b)
+    if !eval_expression
+        eqs, ps, p_ = build_parametrized_eqs(_round!(C*Ξ, digits)', b)
+    else
+        K̃ = _round!(C*Ξ, digits)
+        eqs = K̃*Num[states(b); controls(b)]
+        p_ = []
+        ps = [K̃...]
+    end
 
     # Build the lhs
     if length(eqs) == length(states(b))
@@ -317,18 +327,25 @@ function build_solution(prob::DataDrivenProblem, k, C, B, Q, P, inds, b::Abstrac
         eqs = [d(xs[i]) ~ eq for (i,eq) in enumerate(eqs)]
     end
 
+
     res_ = Koopman(eqs, states(b),
         parameters = [parameters(b); p_],
         controls = controls(b), iv = independent_variable(b),
-        K = k, C = C, Q = Q, P = P, lift = b.f)
+        K = k, C = C, Q = Q, P = P, lift = b.f,
+        eval_expression = eval_expression)
 
-
-    retcode = :success
-    pnew = !isempty(parameters(b)) ? [prob.p; ps] : ps
-    # Equation space
+    retcode = :sucess
     X = prob.DX
-    X_, _, t, U = get_oop_args(prob)
-    Y = res_(X_, pnew, t, U)
+    X_, p_, t, U = get_oop_args(prob)
+
+    pnew = !isempty(parameters(b)) ? [p_; ps] : ps
+    
+    if !eval_expression
+        # Equation space
+        Y = res_(X_, pnew, t, U)
+    else
+        Y = K̃*b(X_, p_, t, U)
+    end
 
     # Build the metrics
     error = norm(X-Y, 2)
@@ -357,5 +374,4 @@ function build_solution(prob::DataDrivenProblem, k, C, B, Q, P, inds, b::Abstrac
     return DataDrivenSolution(
         res_, retcode, pnew, alg, Ξ, inputs, metrics
     )
-    return K
 end
diff --git a/src/solve/koopman.jl b/src/solve/koopman.jl
index 4f532c0f5..8bea60534 100644
--- a/src/solve/koopman.jl
+++ b/src/solve/koopman.jl
@@ -1,6 +1,7 @@
 
 function DiffEqBase.solve(prob::DataDrivenProblem{dType}, alg::AbstractKoopmanAlgorithm;
     B::AbstractArray = [], digits::Int = 10, operator_only::Bool = false,
+    eval_expression = false,
     kwargs...) where {dType <: Number}
     # Check the validity
     @assert is_valid(prob) "The problem seems to be ill-defined. Please check the problem definition."
@@ -39,11 +40,12 @@ function DiffEqBase.solve(prob::DataDrivenProblem{dType}, alg::AbstractKoopmanAl
 
     operator_only && return (K = k, C = C, B = B, Q = Q, P = P)
 
-    return build_solution(prob, k, C, B, Q, P, inds, b, alg, digits = digits)
+    return build_solution(prob, k, C, B, Q, P, inds, b, alg, digits = digits, eval_expression = eval_expression)
 end
 
 function DiffEqBase.solve(prob::DataDrivenProblem{dType}, b::Basis, alg::AbstractKoopmanAlgorithm;
     digits::Int = 10, operator_only::Bool = false,
+    eval_expression = false, 
     kwargs...) where {dType <: Number}
     # Check the validity
     @assert is_valid(prob) "The problem seems to be ill-defined. Please check the problem definition."
@@ -98,5 +100,5 @@ function DiffEqBase.solve(prob::DataDrivenProblem{dType}, b::Basis, alg::Abstrac
 
     operator_only && return (K = k, C = C, B = B, Q = Q, P = P)
 
-    return build_solution(prob, k, C, B, Q, P, inds, b, alg, digits = digits)
+    return build_solution(prob, k, C, B, Q, P, inds, b, alg, digits = digits, eval_expression = eval_expression)
 end
diff --git a/src/solve/sindy.jl b/src/solve/sindy.jl
index b998c0537..0a9bf4e11 100644
--- a/src/solve/sindy.jl
+++ b/src/solve/sindy.jl
@@ -22,7 +22,8 @@ end
 # Main
 function DiffEqBase.solve(p::DataDrivenProblem{dType}, b::Basis, opt::Optimize.AbstractOptimizer;
     normalize::Bool = false, denoise::Bool = false, maxiter::Int = 0,
-    round::Bool = true, kwargs...) where {dType <: Number}
+    round::Bool = true,
+    eval_expression = false, kwargs...) where {dType <: Number}
     # Check the validity
     @assert is_valid(p) "The problem seems to be ill-defined. Please check the problem definition."
 
@@ -48,7 +49,7 @@ function DiffEqBase.solve(p::DataDrivenProblem{dType}, b::Basis, opt::Optimize.A
 
     # Build solution Basis
     return build_solution(
-        p, Ξ, opt, b
+        p, Ξ, opt, b, eval_expression = eval_expression
     )
 end
 
@@ -78,6 +79,7 @@ end
 @views function DiffEqBase.solve(p::DataDrivenProblem{dType}, b::Basis,
     opt::Optimize.AbstractSubspaceOptimizer, implicits::Vector{Num} = Num[];
     normalize::Bool = false, denoise::Bool = false, maxiter::Int = 0,
+    eval_expression = false,
     round::Bool = true, kwargs...) where {dType <: Number}
     # Check the validity
     @assert is_valid(p) "The problem seems to be ill-defined. Please check the problem definition."
@@ -119,6 +121,6 @@ end
     normalize ? rescale_xi!(Ξ, scales, round) : nothing
     # Build solution Basis
     return build_solution(
-        p, Ξ, opt, b, implicits
+        p, Ξ, opt, b, implicits, eval_expression = eval_expression
     )
 end
diff --git a/test/dmd/linear_autonomous.jl b/test/dmd/linear_autonomous.jl
index 568668b22..1e7345788 100644
--- a/test/dmd/linear_autonomous.jl
+++ b/test/dmd/linear_autonomous.jl
@@ -34,7 +34,7 @@ end
     u0 = [10.0; -20.0]
     prob = ODEProblem(f, u0, (0.0, 10.0))
     sol = solve(prob, Tsit5(), saveat = 0.001)
-    
+
     prob = DataDrivenProblem(sol)
 
     for alg in [DMDPINV(), DMDSVD(), TOTALDMD()]
@@ -45,7 +45,7 @@ end
         @test isempty(estimator.B)
         res = solve(prob, alg , operator_only = false)
         m = metrics(res)
-        @test m.Error ./ size(X, 2) < 3e-1
+        @test m.Error ./ size(sol, 2) < 3e-1
         @test Matrix(result(res)) ≈ Matrix(estimator.K)
     end
 end
@@ -74,6 +74,16 @@ end
     end
 end
 
+@testset "Big System" begin
+    # Creates a big system which would resulting in a segfault otherwise
+    X = rand([0, 1], 128, 936);
+    T = collect(LinRange(0, 4.367058580858928, 936));
+    problem = DiscreteDataDrivenProblem(X, T);
+    res1 = solve(problem, DMDSVD(), eval_expression = true)
+    res2 = solve(problem, DMDSVD(), operator_only = true)
+    @test Matrix(result(res1)) == real.(Matrix(res2.K))
+end
+
 # TODO Include the Big System test
 # This fails to generate an equation right now...
 #@testset "Big System" begin
diff --git a/test/dmd/nonlinear_autonomous.jl b/test/dmd/nonlinear_autonomous.jl
index 286d7b8a1..9d02d7aa0 100644
--- a/test/dmd/nonlinear_autonomous.jl
+++ b/test/dmd/nonlinear_autonomous.jl
@@ -21,10 +21,10 @@
     b = result(res)
     m = metrics(res)
     @test isapprox(eigvals(b), [2*p[1]; p[1]; p[2]], atol = 1e-1)
-    @test m.Error/size(X, 2) < 1e-1
+    @test m.Error/size(solution, 2) < 1e-1
 
     _prob = ODEProblem((args...)->b(args...), u0, tspan, parameters(res))
     _sol = solve(_prob, Tsit5(), saveat = solution.t)
-    @test norm(solution - _sol)/size(X, 2) < 1e-1
+    @test norm(solution - _sol)/size(solution, 2) < 1e-1
   end
 end
diff --git a/test/dmd/nonlinear_forced.jl b/test/dmd/nonlinear_forced.jl
index 6ea158e9e..de92689fa 100644
--- a/test/dmd/nonlinear_forced.jl
+++ b/test/dmd/nonlinear_forced.jl
@@ -11,7 +11,7 @@
 
   problem = ODEProblem(slow_manifold, u0, tspan, p)
   solution = solve(problem, Tsit5(), saveat = 0.01)
-  
+
   ufun(u,p,t) = sin(t^2)
 
   prob = ContinuousDataDrivenProblem(solution, U = ufun)
@@ -24,7 +24,7 @@
     b = result(res)
     m = metrics(res)
     @test isapprox(eigvals(b), [2*p[1]; p[1]; p[2]], atol = 1e-1)
-    @test m.Error/size(X, 2) < 3e-1
+    @test m.Error/size(solution, 2) < 3e-1
 
     # TODO This does not work right now, but it should
     #sdict = Dict([y[1] => sin(t^2)])