diff --git a/CondaPkg.toml b/CondaPkg.toml
index 3786d71..39ba32c 100644
--- a/CondaPkg.toml
+++ b/CondaPkg.toml
@@ -1,8 +1,9 @@
 [deps]
-pybamm = "23.1"
 networkx = ""
+casadi = ""
 python = "=3.9"
 
 [pip.deps]
+pybamm = "@ git+https://github.com/pybamm-team/PyBaMM@develop"
 qtconsole = ""
 dfo-ls = ""
diff --git a/pysrc/pack.py b/pysrc/pack.py
index b257157..f39611b 100644
--- a/pysrc/pack.py
+++ b/pysrc/pack.py
@@ -181,7 +181,7 @@ def __init__(
                 sim.built_model.concatenated_rhs, sim.built_model.concatenated_algebraic
             )
 
-        self.timescale = sim.built_model.timescale
+        self.timescale = 1
 
         self.len_cell_algebraic = sim.built_model.len_alg
 
diff --git a/scripts/jl/pack_from_circuit.jl b/scripts/jl/pack_from_circuit.jl
index cc3f8d8..88cd816 100644
--- a/scripts/jl/pack_from_circuit.jl
+++ b/scripts/jl/pack_from_circuit.jl
@@ -41,7 +41,7 @@ else
 end
 
 
-timescale = pyconvert(Float64,pybamm_pack.timescale.evaluate())
+timescale = 1
 cellconverter = pybamm2julia.JuliaConverter(cache_type = "symbolic", inplace=true)
 cellconverter.convert_tree_to_intermediate(pybamm_pack.cell_model)
 cell_str = cellconverter.build_julia_code()
@@ -105,9 +105,22 @@ cells = repeat(vcat(cell_rhs,cell_algebraic),pyconvert(Int, pybamm_pack.num_cell
 thermals = trues(pyconvert(Int,pybamm_pack.len_thermal_eqs))
 differential_vars = vcat(pack_eqs,cells, thermals)
 mass_matrix = sparse(diagm(differential_vars))
-func = ODEFunction(jl_func, mass_matrix=mass_matrix,jac_prototype=jac_sparsity)
+func = ODEFunction(jl_func, mass_matrix=mass_matrix, jac_prototype=jac_sparsity)
 prob = ODEProblem(func, jl_vec, (0.0, 3600/timescale), p)
 
+using IncompleteLU
+function incompletelu(W,du,u,p,t,newW,Plprev,Prprev,solverdata)
+  if newW === nothing || newW
+    Pl = ilu(convert(AbstractMatrix,W), τ = 50.0)
+  else
+    Pl = Plprev
+  end
+  Pl,nothing
+end
+
+
+Base.eltype(::IncompleteLU.ILUFactorization{Tv,Ti}) where {Tv,Ti} = Tv
+
 
-sol = solve(prob, QNDF(linsolve=KLUFactorization(),concrete_jac=true), saveat = collect(range(0,stop=3600/timescale, length=100)))
+sol = solve(prob, QNDF(linsolve=KrylovJL_GMRES(),precs=incompletelu,concrete_jac=true), save_everystep=false)
 
diff --git a/scripts/pack_prototyping.jl b/scripts/pack_prototyping.jl
index b8e6b53..4e3474a 100644
--- a/scripts/pack_prototyping.jl
+++ b/scripts/pack_prototyping.jl
@@ -5,8 +5,8 @@ pack = PyBaMM.pack
 pybamm2julia = PyBaMM.pybamm2julia
 setup_circuit = PyBaMM.setup_circuit
 
-Np = 5
-Ns = 5
+Np = 3
+Ns = 3
 curr = 1.8
 p = nothing 
 t = 0.0
@@ -16,11 +16,11 @@ model = pybamm.lithium_ion.DFN(name="DFN", options=options)
 netlist = setup_circuit.setup_circuit(Np, Ns, I=curr)
 
 
-pybamm_pack = pack.Pack(model, netlist, functional=functional, thermal=true, thermal_type="legacy", top_bc="symmetry")
+pybamm_pack = pack.Pack(model, netlist, functional=functional, thermals=true, thermal_type="legacy", top_bc="symmetry")
 pybamm_pack.build_pack()
 
-#=
-timescale = pyconvert(Float64,pybamm_pack.timescale.evaluate())
+
+timescale = 1
 cellconverter = pybamm2julia.JuliaConverter(cache_type = "symbolic", inplace=true)
 cellconverter.convert_tree_to_intermediate(pybamm_pack.cell_model)
 cell_str = cellconverter.build_julia_code()
@@ -76,4 +76,4 @@ prob = ODEProblem(func, jl_vec, (0.0, 3600/timescale), nothing)
 
 
 sol = solve(prob, QNDF(linsolve=KLUFactorization(),concrete_jac=true), saveat = collect(range(0,stop=3600/timescale, length=100)))
-=#
+
diff --git a/src/diffeq_problems.jl b/src/diffeq_problems.jl
index c861a3d..c9e9aa1 100644
--- a/src/diffeq_problems.jl
+++ b/src/diffeq_problems.jl
@@ -82,7 +82,7 @@ function _problem_setup(sim, tend, inputs;dae_type="implicit",preallocate=true,c
     end
 
     # Scale the time
-    tau = pyconvert(Float64,sim.built_model.timescale.evaluate())
+    tau = 1
     tspan = (0, tend/tau)
     
 
diff --git a/test/test_caches.jl b/test/test_caches.jl
index f1e3ed1..5c57db0 100644
--- a/test/test_caches.jl
+++ b/test/test_caches.jl
@@ -6,6 +6,7 @@ using OrdinaryDiffEq
 
 
 pybamm=pyimport("pybamm")
+np = pyimport("numpy")
 
 @testset "Dual Cache ODE" begin
     # load model
@@ -18,10 +19,10 @@ pybamm=pyimport("pybamm")
 
     # Calculate voltage in Julia
     V = get_variable(sim, sol, "Terminal voltage [V]")
-    t = get_variable(sim, sol, "Time [s]")
+    t = sol.t
 
     # Solve in python
-    sol_pybamm = sim.solve(t)
+    sol_pybamm = sim.solve(np.array(sol.t))
     V_pybamm = pyconvert(Array{Float64},get(sol_pybamm, "Terminal voltage [V]",nothing).data)
     @test all(isapprox.(V_pybamm, V, atol=1e-2))
 end
@@ -39,10 +40,10 @@ end
 
     # Calculate voltage in Julia
     V = get_variable(sim, sol, "Terminal voltage [V]")
-    t = get_variable(sim, sol, "Time [s]")
+    t = sol.t
 
     # Solve in python
-    sol_pybamm = sim.solve(t)
+    sol_pybamm = sim.solve(np.array(sol.t))
     V_pybamm = pyconvert(Array{Float64},get(sol_pybamm, "Terminal voltage [V]",nothing).data)
     @test all(isapprox.(V_pybamm, V, atol=1e-2))
 end
@@ -59,10 +60,10 @@ end
 
     # Calculate voltage in Julia
     V = get_variable(sim, sol, "Terminal voltage [V]")
-    t = get_variable(sim, sol, "Time [s]")
+    t = sol.t
 
     # Solve in python
-    sol_pybamm = sim.solve(t)
+    sol_pybamm = sim.solve(np.array(sol.t))
     V_pybamm = pyconvert(Array{Float64},get(sol_pybamm, "Terminal voltage [V]",nothing).data)
     @test all(isapprox.(V_pybamm, V, atol=1e-2))
 end
@@ -78,10 +79,10 @@ end
 
     # Calculate voltage in Julia
     V = get_variable(sim, sol, "Terminal voltage [V]")
-    t = get_variable(sim, sol, "Time [s]")
+    t = sol.t
 
     # Solve in python
-    sol_pybamm = sim.solve(t)
+    sol_pybamm = sim.solve(np.array(sol.t))
     V_pybamm = pyconvert(Array{Float64},get(sol_pybamm, "Terminal voltage [V]",nothing).data)
     @test all(isapprox.(V_pybamm, V, atol=1e-2))
 end
diff --git a/test/test_events.jl b/test/test_events.jl
index 0a5bbc2..3b76632 100644
--- a/test/test_events.jl
+++ b/test/test_events.jl
@@ -9,7 +9,7 @@ pybamm = pyimport("pybamm")
     model = pybamm.lithium_ion.SPMe(name="SPMe")
     sim = pybamm.Simulation(model)
     prob,cbs = get_ode_problem(sim)
-    prob = remake(prob,tspan=(0,0.5))
+    prob = remake(prob,tspan=(0,7200))
     event_to_test = sim.built_model.events[3]
     problem_size = length(prob.u0)
     sol = solve(prob,Rodas5(autodiff=false),callback=cbs)
diff --git a/test/test_full_models.jl b/test/test_full_models.jl
index f01c03d..b9f4fe0 100644
--- a/test/test_full_models.jl
+++ b/test/test_full_models.jl
@@ -5,6 +5,7 @@ using Sundials
 using OrdinaryDiffEq
 
 pybamm = pyimport("pybamm")
+np = pyimport("numpy")
 
 @testset "Compare with PyBaMM: SPM" begin
     # load model
@@ -17,7 +18,7 @@ pybamm = pyimport("pybamm")
 
     # Calculate voltage in Julia
     V = get_variable(sim, sol, "Terminal voltage [V]")
-    t = get_variable(sim, sol, "Time [s]")
+    t = np.array(sol.t)
 
     # Solve in python
     sol_pybamm = sim.solve(t)
@@ -36,12 +37,12 @@ end
 
     # Calculate voltage in Julia
     V = get_variable(sim, sol, "Terminal voltage [V]")
-    t = get_variable(sim, sol, "Time [s]")
+    t = np.array(sol.t)
 
     # Solve in python
     sol_pybamm = sim.solve(t)
     V_pybamm = pyconvert(Array{Float64},get(sol_pybamm, "Terminal voltage [V]",nothing).data)
-    @test all(isapprox.(V_pybamm, V, atol=1e-4))
+    @test all(isapprox.(V_pybamm, V, atol=1e-3))
 end
 
 @testset "Compare with PyBaMM: DFN" begin
@@ -55,7 +56,7 @@ end
 
     # Calculate voltage in Julia
     V = get_variable(sim, sol, "Terminal voltage [V]")
-    t = get_variable(sim, sol, "Time [s]")
+    t = np.array(sol.t)
 
     # Solve in python
     sol_pybamm = sim.solve(t)
@@ -75,7 +76,7 @@ end
 
     # Calculate voltage in Julia
     V = get_variable(sim, sol, "Terminal voltage [V]")
-    t = get_variable(sim, sol, "Time [s]")
+    t = np.array(sol.t)
 
     # Solve in python
     sol_pybamm = sim.solve(t)
diff --git a/test/test_pack.jl b/test/test_pack.jl
index a276428..4d92608 100644
--- a/test/test_pack.jl
+++ b/test/test_pack.jl
@@ -30,7 +30,7 @@ using PyBaMM
   pybamm_pack.build_pack()
 
 
-  timescale = pyconvert(Float64,pybamm_pack.timescale.evaluate())
+  timescale = 1
   cellconverter = pybamm2julia.JuliaConverter(cache_type = "symbolic", inplace=true)
   cellconverter.convert_tree_to_intermediate(pybamm_pack.cell_model)
   cell_str = cellconverter.build_julia_code()
@@ -122,7 +122,7 @@ end
 
   pybamm_pack = pack.Pack(model, circuit_graph, functional=functional, thermals=thermals, operating_mode = "CV")
   pybamm_pack.build_pack()
-  timescale = pyconvert(Float64,pybamm_pack.timescale.evaluate())
+  timescale = 1
   cellconverter = pybamm2julia.JuliaConverter(cache_type = "symbolic", inplace=true)
   cellconverter.convert_tree_to_intermediate(pybamm_pack.cell_model)
   cell_str = cellconverter.build_julia_code()