diff --git a/Project.toml b/Project.toml
index 80681a8..1cb8a09 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,15 +1,15 @@
 name = "PETLION"
 uuid = "5e0a28e4-193c-47fa-bbb8-9c901cc1ac2c"
 authors = ["Marc D. Berliner", "Richard D. Braatz"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 GeneralizedGenerated = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
@@ -18,6 +18,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
@@ -35,7 +36,7 @@ SpecialFunctions = "1"
 StatsBase = "0.3 - 0.33"
 Sundials = "4"
 Symbolics = "0.1, 1, 2, 3"
-julia = "1, 1.5, 1.6, 1.7"
+julia = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/src/PETLION.jl b/src/PETLION.jl
index d2e0a1d..9aacbed 100644
--- a/src/PETLION.jl
+++ b/src/PETLION.jl
@@ -14,6 +14,7 @@ using RecipesBase
 using SpecialFunctions: erf
 using SHA: sha1
 
+import SuiteSparse
 import Sundials
 import LinearAlgebra
 
diff --git a/src/checks.jl b/src/checks.jl
index 0cf3697..f5f731f 100644
--- a/src/checks.jl
+++ b/src/checks.jl
@@ -1,6 +1,6 @@
-@inline function check_simulation_stop!(model, t::Float64, Y, YP, run::AbstractRun, p, bounds, opts::options_model;
+@inline function check_simulation_stop!(model, t::Float64, Y, YP, run::AbstractRun, p, bounds, opts::options_model_immutable{T};
     ϵ::Float64 = t < 1.0 ? opts.reltol : 0.0,
-    )
+    ) where T<:Function
     
     if t ≥ run.tf
         run.info.flag = 0
@@ -21,6 +21,7 @@
     check_stop_c_e(       p, run, model, Y, YP, bounds, ϵ, I)
     check_stop_η_plating( p, run, model, Y, YP, bounds, ϵ, I)
     check_stop_dfilm(     p, run, model, Y, YP, bounds, ϵ, I)
+    opts.stop_function(   p, run, model, Y, YP, bounds, ϵ, I)
 
     return nothing
 end
@@ -220,7 +221,7 @@ end
 
 
 
-@inline function check_solve(run::run_constant, model::R1, int::R2, p, bounds, opts::R5, funcs, keep_Y::Bool, iter::Int64, Y::Vector{Float64}, t::Float64) where {R1<:model_output,R2<:Sundials.IDAIntegrator,R5<:options_model}
+@inline function check_solve(run::run_constant, model::R1, int::R2, p, bounds, opts::R5, funcs, keep_Y::Bool, iter::Int64, Y::Vector{Float64}, t::Float64) where {R1<:model_output,R2<:Sundials.IDAIntegrator,R5<:AbstractOptionsModel}
     if t === int.tprev
         # Sometimes the initial step at t = 0 can be too large. This reduces the step size
         if t === 0.0
@@ -245,7 +246,7 @@ end
     return true
 end
 
-@inline function check_solve(run::run_function, model::R1, int::R2, p::param, bounds::boundary_stop_conditions, opts::R5, funcs, keep_Y::Bool, iter::Int64, Y::Vector{Float64}, t::Float64) where {R1<:model_output,R2<:Sundials.IDAIntegrator,R5<:options_model}
+@inline function check_solve(run::run_function, model::R1, int::R2, p::param, bounds::boundary_stop_conditions, opts::R5, funcs, keep_Y::Bool, iter::Int64, Y::Vector{Float64}, t::Float64) where {R1<:model_output,R2<:Sundials.IDAIntegrator,R5<:AbstractOptionsModel}
     if iter === opts.maxiters
         error("Reached max iterations of $(opts.maxiters) at t = $(int.t)")
     elseif within_bounds(run)
@@ -264,7 +265,7 @@ end
     end
 end
 
-@inline function check_reinitialization!(model::R1, int::R2, run::R3, p::R4, bounds::R5, opts::R6, funcs) where {R1<:model_output, R2<:Sundials.IDAIntegrator, R3<:AbstractRun,R4<:param,R5<:boundary_stop_conditions,R6<:options_model}
+@inline function check_reinitialization!(model::R1, int::R2, run::R3, p::R4, bounds::R5, opts::R6, funcs) where {R1<:model_output, R2<:Sundials.IDAIntegrator, R3<:AbstractRun,R4<:param,R5<:boundary_stop_conditions,R6<:AbstractOptionsModel}
     """
     Checking the current function for discontinuities.
     If there is a significant change in current after a step size of dt = reltol,
@@ -273,21 +274,23 @@ end
     
     Y = int.u.v
     YP = int.du.v
+    # take a step of Δt = the relative tolerance
     t_new = int.t + opts.reltol
 
     value_old = value(run)
     value_new = run.func(t_new,Y,YP,p)
     
+    # if the function values at t vs. t + Δt are very different (i.e., there is a discontinuity)
+    # then reinitialize the DAE at t + Δt
     if !≈(value_old, value_new, atol=opts.abstol, rtol=opts.reltol)
         initialize_states!(p,Y,YP,run,opts,funcs,(@inbounds model.SOC[end]); t=t_new)
-        #run.value .= value_new
 
         Sundials.IDAReInit(int.mem, t_new, Y, YP)
     end
     return nothing
 end
 
-@inline function check_errors_parameters_runtime(p::R1,opts::R2) where {R1<:param,R2<:options_model}
+@inline function check_errors_parameters_runtime(p::R1,opts::R2) where {R1<:param,R2<:AbstractOptionsModel}
     ϵ_sp, ϵ_sn = active_material(p)
 
     if ( ϵ_sp > 1 ) error("ϵ_p + ϵ_fp must be ∈ [0, 1)") end
@@ -307,4 +310,4 @@ function check_errors_initial(θ, numerics, N)
 end
 
 check_is_hold(x::Symbol,model::model_output) = (x===:hold) && (!isempty(model) ? true : error("Cannot use `:hold` without a previous model."))
-check_is_hold(x,model) = false
+check_is_hold(::Any,::model_output) = false
\ No newline at end of file
diff --git a/src/external.jl b/src/external.jl
index f0b33f0..2182c02 100644
--- a/src/external.jl
+++ b/src/external.jl
@@ -413,7 +413,11 @@ end
 
 function strings_directory_func(N::discretizations_per_section, numerics::T; create_dir=false) where T<:options_numerical
 
-    dir_saved_models = joinpath(options[:FILE_DIRECTORY], "saved_models")
+    dir_saved_models = "saved_models"
+    # If the file directory is not specified, use the current working directory
+    if !isnothing(options[:FILE_DIRECTORY])
+        dir_saved_models = joinpath(options[:FILE_DIRECTORY], dir_saved_models)
+    end
     
     if create_dir && !isdir(dir_saved_models)
         mkdir(dir_saved_models)
diff --git a/src/generate_functions.jl b/src/generate_functions.jl
index 9bf6558..79ff3aa 100644
--- a/src/generate_functions.jl
+++ b/src/generate_functions.jl
@@ -1,7 +1,8 @@
-const PETLION_VERSION = (0,2,0)
+const PETLION_VERSION = (0,2,1)
 const options = Dict{Symbol,Any}(
     :SAVE_SYMBOLIC_FUNCTIONS => true,
-    :FILE_DIRECTORY => pwd(),
+    :FILE_DIRECTORY => nothing,
+    :FACTORIZATION_METHOD => :KLU, # :KLU or :LU
 )
 
 function load_functions(p::AbstractParam)
diff --git a/src/model_evaluation.jl b/src/model_evaluation.jl
index fee616f..009e106 100644
--- a/src/model_evaluation.jl
+++ b/src/model_evaluation.jl
@@ -5,21 +5,23 @@
     I   = nothing, # constant current C-rate. also accepts :rest
     V   = nothing, # constant voltage. also accepts :hold if continuing simulation
     P   = nothing, # constant power.   also accepts :hold if continuing simulation
-    SOC::Number   = p.opts.SOC, # initial SOC of the simulation. only valid if not continuing simulation
-    outputs::R2   = p.opts.outputs, # model output states
-    abstol        = p.opts.abstol, # absolute tolerance in DAE solve
-    reltol        = p.opts.reltol, # relative tolerance in DAE solve
-    abstol_init   = abstol, # absolute tolerance in initialization
-    reltol_init   = reltol, # relative tolerance in initialization
-    maxiters      = p.opts.maxiters, # maximum solver iterations
-    check_bounds  = p.opts.check_bounds, # check if the boundaries (V_min, SOC_max, etc.) are satisfied
-    reinit        = p.opts.reinit, # reinitialize the initial guess
-    verbose       = p.opts.verbose, # print information about the run
-    interp_final  = p.opts.interp_final, # interpolate the final points if a boundary is hit
-    tstops        = p.opts.tstops, # times the solver explicitly solves for
-    tdiscon       = p.opts.tdiscon, # times of known discontinuities in the current function
-    interp_bc     = p.opts.interp_bc, # :interpolate or :extrapolate
-    save_start    = p.opts.save_start, # warm-start for the initial guess
+    SOC::Number     = p.opts.SOC, # initial SOC of the simulation. only valid if not continuing simulation
+    outputs::R2     = p.opts.outputs, # model output states
+    abstol          = p.opts.abstol, # absolute tolerance in DAE solve
+    reltol          = p.opts.reltol, # relative tolerance in DAE solve
+    abstol_init     = abstol, # absolute tolerance in initialization
+    reltol_init     = reltol, # relative tolerance in initialization
+    maxiters        = p.opts.maxiters, # maximum solver iterations
+    check_bounds    = p.opts.check_bounds, # check if the boundaries (V_min, SOC_max, etc.) are satisfied
+    reinit          = p.opts.reinit, # reinitialize the initial guess
+    verbose         = p.opts.verbose, # print information about the run
+    interp_final    = p.opts.interp_final, # interpolate the final points if a boundary is hit
+    tstops          = p.opts.tstops, # times the solver explicitly solves for
+    tdiscon         = p.opts.tdiscon, # times of known discontinuities in the current function
+    interp_bc       = p.opts.interp_bc, # :interpolate or :extrapolate
+    save_start      = p.opts.save_start, # warm-start for the initial guess
+    stop_function   = p.opts.stop_function,
+    calc_integrator = p.opts.calc_integrator,
     V_max         = p.bounds.V_max,
     V_min         = p.bounds.V_min,
     SOC_max       = p.bounds.SOC_max,
@@ -59,7 +61,7 @@
     run = get_run(I,V,P,t0,tf,p,model)
     
     # putting opts and bounds into a structure. 
-    opts = options_model(outputs, Float64(SOC), abstol, reltol, abstol_init, reltol_init, maxiters, check_bounds, reinit, verbose, interp_final, tstops, tdiscon, interp_bc, save_start, var_keep)
+    opts = options_model_immutable(outputs, Float64(SOC), abstol, reltol, abstol_init, reltol_init, maxiters, check_bounds, reinit, verbose, interp_final, tstops, tdiscon, interp_bc, save_start, var_keep, stop_function, calc_integrator)
     bounds = boundary_stop_conditions(V_max, V_min, SOC_max, SOC_min, T_max, c_s_n_max, I_max, I_min, η_plating_min, c_e_min, dfilm_max)
     
     # getting the initial conditions and run setup
@@ -146,7 +148,7 @@ end
 
 @inline within_bounds(run::AbstractRun) = run.info.flag === -1
 
-@inline function initialize_model!(model::model_struct, p::param{jac}, run::T, bounds::boundary_stop_conditions, opts::options_model, res_I_guess=1.0) where {jac<:AbstractJacobian,method<:AbstractMethod,T<:AbstractRun{method,<:Any},model_struct<:Union{model_output,Vector{Float64}}}
+@inline function initialize_model!(model::model_struct, p::param{jac}, run::T, bounds::boundary_stop_conditions, opts::AbstractOptionsModel, res_I_guess=nothing) where {jac<:AbstractJacobian,method<:AbstractMethod,T<:AbstractRun{method,<:Any},model_struct<:Union{model_output,Vector{Float64}}}
     if !haskey(p.funcs,run)
         get_method_funcs!(p,run)
         funcs = p.funcs(run)
@@ -198,13 +200,13 @@ end
     return int, funcs, model
 end
 
-@inline function retrieve_integrator(run::T, p::param, funcs::Jac_and_res{<:Sundials.IDAIntegrator}, Y0, YP0, opts::options_model, new_run::Bool) where {method<:AbstractMethod,T<:AbstractRun{method,<:Any}}
+@inline function retrieve_integrator(run::T, p::param, funcs::Jac_and_res{<:Sundials.IDAIntegrator}, Y0, YP0, opts::AbstractOptionsModel, new_run::Bool) where {method<:AbstractMethod,T<:AbstractRun{method,<:Any}}
     """
     If the model has previously been evaluated for a constant run simulation, you can reuse
     the integrator function with its cache instead of creating a new one
     """
     
-    if isempty(funcs.int)
+    if isempty(funcs.int) || opts.calc_integrator
         int = create_integrator(run, p, funcs, Y0, YP0, opts)
     else
         # reuse the integrator cache
@@ -223,7 +225,7 @@ end
     return int
 end
 
-@inline function create_integrator(run::T, p::param, funcs::Q, Y0, YP0, opts::options_model) where {T<:AbstractRun,Q<:Jac_and_res}
+@inline function create_integrator(run::T, p::param, funcs::Q, Y0, YP0, opts::AbstractOptionsModel) where {T<:AbstractRun,Q<:Jac_and_res}
     R_full = funcs.R_full
     J_full = funcs.J_full
     DAEfunc = DAEFunction(
@@ -242,7 +244,7 @@ end
 
     return int
 end
-@inline function postfix_integrator!(int::Sundials.IDAIntegrator, run::AbstractRun, opts::options_model, new_run::Bool)
+@inline function postfix_integrator!(int::Sundials.IDAIntegrator, run::AbstractRun, opts::AbstractOptionsModel, new_run::Bool)
     tstops = int.opts.tstops.valtree
     empty!(tstops)
     
@@ -258,7 +260,7 @@ end
     return nothing
 end
 
-@inline function solve!(model,int::R1,run::R2,p,bounds,opts::R3,funcs) where {R1<:Sundials.IDAIntegrator,R2<:AbstractRun,R3<:options_model}
+@inline function solve!(model,int::R1,run::R2,p,bounds,opts::R3,funcs) where {R1<:Sundials.IDAIntegrator,R2<:AbstractRun,R3<:AbstractOptionsModel}
     keep_Y = opts.var_keep.Y
     Y  = int.u.v
     YP = int.du.v
@@ -281,10 +283,10 @@ end
     return nothing
 end
 
-@inline function exit_simulation!(p::R1, model::R2, run::R3, bounds::R4, int::R5, opts::R6; cancel_interp::Bool=false) where {R1<:param,R2<:model_output,R3<:AbstractRun,R4<:boundary_stop_conditions,R5<:Sundials.IDAIntegrator,R6<:options_model}
+@inline function exit_simulation!(p::R1, model::R2, run::R3, bounds::R4, int::R5, opts::R6; cancel_interp::Bool=false) where {R1<:param,R2<:model_output,R3<:AbstractRun,R4<:boundary_stop_conditions,R5<:Sundials.IDAIntegrator,R6<:AbstractOptionsModel}
     # if a stop condition (besides t = tf) was reached
     if !cancel_interp
-        if opts.interp_final && !(run.info.flag === 0)
+        if opts.interp_final && !(run.info.flag === 0) && int.t > 1
             interp_final_points!(p, model, run, bounds, int, opts)
         else
             set_var!(model.Y, opts.var_keep.Y ? copy(int.u.v) : int.u.v, opts.var_keep.Y)
@@ -315,8 +317,13 @@ end
     return nothing
 end
 
-@views @inbounds @inline function interp_final_points!(p::R1, model::R2, run::R3, bounds::R4, int::R5, opts::R6) where {R1<:param,R2<:model_output,R3<:AbstractRun,R4<:boundary_stop_conditions,R5<:Sundials.IDAIntegrator,R6<:options_model}
-    YP = opts.var_keep.YP ? model.YP[end] : Float64[]
+@views @inbounds @inline function interp_final_points!(p::R1, model::R2, run::R3, bounds::R4, int::R5, opts::R6) where {R1<:param,R2<:model_output,R3<:AbstractRun,R4<:boundary_stop_conditions,R5<:Sundials.IDAIntegrator,R6<:AbstractOptionsModel}
+    if opts.var_keep.YP
+        YP = length(model.YP) > 1 ? bounds.t_final_interp_frac.*(model.YP[end] .- model.YP[end-1]) .+ model.YP[end-1] : model.YP[end]
+    else
+        YP = Float64[]
+    end
+    
     t = bounds.t_final_interp_frac*(int.t - int.tprev) + int.tprev
     
     set_var!(model.Y,  bounds.t_final_interp_frac.*(int.u.v .- model.Y[end]) .+ model.Y[end], opts.var_keep.Y)
@@ -341,7 +348,7 @@ end
     
     return key, key_exists
 end
-@inline function initialize_states!(p::param, Y0::T, YP0::T, run::AbstractRun, opts::options_model, funcs::Jac_and_res, SOC::Float64;kw...) where {T<:Vector{Float64}}
+@inline function initialize_states!(p::param, Y0::T, YP0::T, run::AbstractRun, opts::AbstractOptionsModel, funcs::Jac_and_res, SOC::Float64;kw...) where {T<:Vector{Float64}}
     if opts.save_start
         key, key_exists = save_start_init!(Y0, run, p, SOC)
         
@@ -357,7 +364,10 @@ end
     return nothing
 end
 
-@inline function newtons_method!(p::param,Y::R1,YP::R1,run,opts::options_model,R_alg::T1,R_diff::T2,J_alg::T3;
+
+@inline factorization!(L::SuiteSparse.UMFPACK.UmfpackLU{Float64, Int64},A::SparseMatrixCSC{Float64, Int64}) = LinearAlgebra.lu!(L,A)
+@inline factorization!(L::KLUFactorization{Float64, Int64},A::SparseMatrixCSC{Float64, Int64}) = klu!(L,A)
+@inline function newtons_method!(p::param,Y::R1,YP::R1,run,opts::AbstractOptionsModel,R_alg::T1,R_diff::T2,J_alg::T3;
     itermax::Int64=100, t::Float64=0.0
     ) where {R1<:Vector{Float64},T1<:residual_combined,T2<:residual_combined,T3<:jacobian_combined}
 
@@ -367,14 +377,14 @@ end
     YP   .= 0.0
     J     = J_alg.sp
     γ     = 0.0
-    L     = J_alg.L # KLU factorization
+    L     = J_alg.L # factorization
     
     # starting loop for Newton's method
     @inbounds for iter in 1:itermax
         # updating res, Y, and J
         R_alg(res,t,Y,YP,p,run)
         J_alg(t,Y,YP,γ,p,run)
-        klu!(L, J)
+        factorization!(L, J)
         
         Y_old .= Y_new
         Y_new .-= L\res
diff --git a/src/outputs.jl b/src/outputs.jl
index e31e856..eaedc8c 100644
--- a/src/outputs.jl
+++ b/src/outputs.jl
@@ -1,20 +1,21 @@
+const empty_vector_of_array = VectorOfArray(Array{Vector{Float64}}(Float64[]))
 Base.@kwdef struct model_states{vec1D,vec2D,R1}
     """
     If you want to add anything to this struct, you must also check/modify set_vars!`.
     Otherwise, it may not work as intended
     """
     # Matrices (vectors in space and time)
-    Y::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
-    YP::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
-    c_e::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
-    c_s_avg::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
-    T::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
-    film::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
-    Q::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
-    j::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
-    j_s::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
-    Φ_e::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
-    Φ_s::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
+    Y::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
+    YP::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
+    c_e::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
+    c_s_avg::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
+    T::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
+    film::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
+    Q::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
+    j::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
+    j_s::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
+    Φ_e::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
+    Φ_s::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
     # Vectors (vectors in time, not space)
     I::vec1D = ( vec1D === Array{Float64,1} ) ? Float64[] : nothing
     t::vec1D = ( vec1D === Array{Float64,1} ) ? Float64[] : nothing
diff --git a/src/physics_equations/auxiliary_states_and_coefficients.jl b/src/physics_equations/auxiliary_states_and_coefficients.jl
index 3d210a7..ad57c65 100644
--- a/src/physics_equations/auxiliary_states_and_coefficients.jl
+++ b/src/physics_equations/auxiliary_states_and_coefficients.jl
@@ -489,10 +489,13 @@ function limiting_electrode(p::AbstractParam)
     θ = p.θ
     ϵ_sp, ϵ_sn = active_material(p)
 
-    if ϵ_sp*θ[:l_p]*θ[:c_max_p]*(θ[:θ_min_p] - θ[:θ_max_p]) > ϵ_sn*θ[:l_n]*θ[:c_max_n]*(θ[:θ_max_n] - θ[:θ_min_n])
-        return :p
+    Q_p = ϵ_sp*θ[:l_p]*θ[:c_max_p]*(θ[:θ_min_p] - θ[:θ_max_p])
+    Q_n = ϵ_sn*θ[:l_n]*θ[:c_max_n]*(θ[:θ_max_n] - θ[:θ_min_n])
+
+    if Q_p > Q_n
+        return :p, Q_p
     else
-        return :n
+        return :n, Q_n
     end
 end
 
@@ -501,8 +504,8 @@ end
     """
     Calculate the 1C current density (A⋅hr/m²) based on the limiting electrode
     """
-    F = 96485.3321233
-
+    F = const_Faradays
+    
     ϵ_sp = 1.0 - (θ[:ϵ_fp] + θ[:ϵ_p])
     ϵ_sn = 1.0 - (θ[:ϵ_fn] + θ[:ϵ_n])
 
diff --git a/src/physics_equations/scalar_residual.jl b/src/physics_equations/scalar_residual.jl
index dd41daf..4df893c 100644
--- a/src/physics_equations/scalar_residual.jl
+++ b/src/physics_equations/scalar_residual.jl
@@ -314,13 +314,25 @@ function combine_Jac_and_res(p,J_sp_base,J_base_func,J_sp_scalar,J_scalar_func,
     return Jac_and_res(J_full,R_full,J_alg,R_diff,R_alg,Sundials.IDAIntegrator[])
 end
 
+function factorization(x...;kw...)
+    method =  Symbol(lowercase(String(options[:FACTORIZATION_METHOD])))
+    if method === :lu
+        L = LinearAlgebra.lu(x...;kw...)
+    elseif method === :klu
+        L = klu(x...;kw...)
+    else
+        error("FACTORIZATION_METHOD $(options[:FACTORIZATION_METHOD]) is not supported.")
+    end
+    return L
+end
 function _get_jacobian_combined(J_sp_base,J_base_func,J_sp_scalar,J_scalar_func,θ_tot,θ_keys;lu_decomposition=false)
     J_sp = [J_sp_base; J_sp_scalar']
 
     if lu_decomposition
-        L = klu(J_sp)
+        J_sp.nzval .= rand(length(J_sp.nzval))
+        L = factorization(J_sp)
     else
-        L = klu(sparse([1],[1],[1.0]))
+        L = factorization(sparse([1],[1],[1.0]))
     end
 
     ind_base   = findall(J_sp.rowval .< size(J_sp,1))
diff --git a/src/set_variables.jl b/src/set_variables.jl
index 8134124..734f103 100644
--- a/src/set_variables.jl
+++ b/src/set_variables.jl
@@ -8,7 +8,7 @@
         R3<:Vector{Float64},
         R4<:Float64,
         R5<:AbstractRun,
-        R6<:options_model,
+        R6<:AbstractOptionsModel,
         R7<:boundary_stop_conditions,
         R8<:Function,
         }
diff --git a/src/structures.jl b/src/structures.jl
index 2b7e9a7..e6c32c1 100644
--- a/src/structures.jl
+++ b/src/structures.jl
@@ -107,6 +107,7 @@ struct jacobian_combined{
     T1<:Function,
     T2<:Union{SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false},SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}},
     T3<:Function,
+    T4<:LinearAlgebra.Factorization{Float64},
     }
     sp::SparseMatrixCSC{Float64,Int64}
     base_func::T1
@@ -115,7 +116,7 @@ struct jacobian_combined{
     J_scalar::SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false}
     θ_tot::Vector{Float64}
     θ_keys::Vector{Symbol}
-    L::KLUFactorization{Float64, Int64}
+    L::T4
 end
 Base.getindex(J::jacobian_combined,i...) = getindex(J.sp,i...)
 Base.axes(J::jacobian_combined,i...) = axes(J.sp,i...)
@@ -226,7 +227,8 @@ const indices_states = model_states{
     Nothing,
 }
 
-Base.@kwdef mutable struct options_model
+abstract type AbstractOptionsModel end
+Base.@kwdef mutable struct options_model <: AbstractOptionsModel
     outputs::Tuple = (:t, :V)
     SOC::Number = 1.0
     abstol::Float64 = 1e-6
@@ -243,7 +245,31 @@ Base.@kwdef mutable struct options_model
     interp_bc::Symbol = :interpolate
     save_start::Bool = false
     var_keep::states_logic = model_states_logic()
-end
+    stop_function::Function = (x...) -> nothing
+    calc_integrator::Bool = false
+end
+struct options_model_immutable{T<:Function} <: AbstractOptionsModel
+    outputs::Tuple
+    SOC::Number
+    abstol::Float64
+    reltol::Float64
+    abstol_init::Float64
+    reltol_init::Float64
+    maxiters::Int64
+    check_bounds::Bool
+    reinit::Bool
+    verbose::Bool
+    interp_final::Bool
+    tstops::Vector{Float64}
+    tdiscon::Vector{Float64}
+    interp_bc::Symbol
+    save_start::Bool
+    var_keep::states_logic
+    stop_function::T
+    calc_integrator::Bool
+end
+# convert from the mutable to immutable struct
+options_model_immutable(opts::options_model) = options_model_immutable((getproperty(opts,field) for field in fieldnames(options_model))...)
 
 struct save_start_info{T<:AbstractMethod}
     method::T
@@ -282,7 +308,6 @@ model_funcs(x...) = model_funcs(x...,
     Dict{DataType,Dict{DataType,Jac_and_res}}(),
     Dict{DataType,Jac_and_res}()
     )
-Base.empty!(f::model_funcs) = ([empty!(getproperty(f,field)) for (_type,field) in zip(fieldtypes(model_funcs),fieldnames(model_funcs)) if _type <: Dict];nothing)
 
 struct param{T<:AbstractJacobian,temp,solid_diff,Fickian,age} <: AbstractParam{T,temp,solid_diff,Fickian,age}
     θ::Dict{Symbol,Float64}
@@ -323,7 +348,7 @@ struct run_results{T<:AbstractRun}
     info::run_info
     run_index::UnitRange{Int64}
     int::Sundials.IDAIntegrator
-    opts::options_model
+    opts::options_model_immutable
     bounds::boundary_stop_conditions
     N::discretizations_per_section
     numerics::options_numerical
@@ -349,6 +374,9 @@ end
 Base.lastindex(model::T) where T<:model_output = length(model)
 Base.firstindex(::T) where T<:model_output = 1
 
+Base.empty!(f::model_funcs) = ([empty!(getproperty(f,field)) for (_type,field) in zip(fieldtypes(model_funcs),fieldnames(model_funcs)) if _type <: Dict];nothing)
+Base.empty!(p::param) = empty!(p.funcs)
+
 
 ## Modifying Base functions
 @recipe function plot(model::model_output, x_name::Symbol=:V;linewidth=2,legend=false)
@@ -438,7 +466,6 @@ function C_rate_string(I::Number;digits::Int64=4)
 end
 
 function Base.show(io::IO, p::AbstractParam)
-    
     function show_bounds(title, min, max, units="")
         if isnan(min) && isnan(max)
             return ""