diff --git a/Project.toml b/Project.toml
index f93d858..0cebd15 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ExtendableASGFEM"
 uuid = "a4750d0b-645d-4217-869c-e25952806e24"
-version = "1.0.0"
+version = "1.0.1"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Martin Eigel <martin.eigel@wias-berlin.de>"]
 
 [deps]
@@ -16,6 +16,7 @@ IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -30,15 +31,15 @@ Distributions = "0.25.120"
 DocStringExtensions = "0.9.4"
 DoubleFloats = "1.4.3"
 ExplicitImports = "1"
-ExtendableFEM = "1.4"
-ExtendableFEMBase = "1.4"
+ExtendableFEM = "1.10"
+ExtendableFEMBase = "1.5.1"
 ExtendableGrids = "1.13"
-ExtendableSparse = "1.7"
+ExtendableSparse = "2.0.1"
 GridVisualize = "1.12"
 IterativeSolvers = "0.9.2"
 Krylov = "0.10.1"
-LaTeXStrings = "1.4.0"
 LinearAlgebra = "1.9"
+LinearSolve = "3.66.0"
 Printf = "1.9"
 Random = "1.9"
 SparseArrays = "1.9"
diff --git a/src/ExtendableASGFEM.jl b/src/ExtendableASGFEM.jl
index 382004d..c940957 100644
--- a/src/ExtendableASGFEM.jl
+++ b/src/ExtendableASGFEM.jl
@@ -24,13 +24,13 @@ using ExtendableGrids: ExtendableGrids, Adjacency, BFaceFaces, CellFaces,
     atranspose, eval_trafo!, interpolate!,
     max_num_targets_per_source, num_cells, num_nodes,
     num_sources, num_targets, unique, update_trafo!
-using ExtendableSparse: ExtendableSparse, ExtendableSparseMatrix,
-    LUFactorization, flush!
+using ExtendableSparse: ExtendableSparse, ExtendableSparseMatrix, flush!
 using GridVisualize: GridVisualize, GridVisualizer, scalarplot, scalarplot!
 using IterativeSolvers: IterativeSolvers
 using Krylov: Krylov
 using LinearAlgebra: LinearAlgebra, SymTridiagonal, dot, eigvals, eigvecs,
-    ldiv!, mul!, norm
+    ldiv!, mul!, norm, lu
+using LinearSolve: LinearSolve, LUFactorization
 using Printf: Printf, @printf
 using Random: Random, rand!
 using SparseArrays: SparseArrays, SparseMatrixCSC, nzrange, rowvals
diff --git a/src/modelproblems/solvers_logpoisson_dual.jl b/src/modelproblems/solvers_logpoisson_dual.jl
index 5651a04..026e06a 100644
--- a/src/modelproblems/solvers_logpoisson_dual.jl
+++ b/src/modelproblems/solvers_logpoisson_dual.jl
@@ -22,7 +22,7 @@ struct MySystemLogDual{Tv, MT, VT, GT}
 end
 Base.size(S::MySystemLogDual) = S.nmodes .* (size(S.A.entries) .+ size(S.B.entries)[2])
 
-struct MyPreconditionerLogDual{Tv, FAC <: ExtendableSparse.AbstractFactorization}
+struct MyPreconditionerLogDual{Tv, FAC}
     LUS::FAC
     DA::Array{Tv, 1}
     temp::Array{Tv, 1}
@@ -55,7 +55,7 @@ function MyPreconditionerLogDual(A::ExtendableSparseMatrix{Tv, Ti}, B::Extendabl
 
     # compute LU factorisation of S
     flush!(S)
-    LUS = LUFactorization(S)
+    LUS = lu(S.cscmatrix)
 
     # temporary storage array for solver
     temp = zeros(Tv, size(B, 2))
diff --git a/src/modelproblems/solvers_logpoisson_primal.jl b/src/modelproblems/solvers_logpoisson_primal.jl
index d6786f9..c78d794 100644
--- a/src/modelproblems/solvers_logpoisson_primal.jl
+++ b/src/modelproblems/solvers_logpoisson_primal.jl
@@ -34,18 +34,16 @@ function MyPreconditionerLogPrimal(A::ExtendableSparseMatrix{Tv, Ti}, bdofs, nmo
         DA[j] = A[j, j]
     end
 
-    # compute LU factorisation of S
+    # compute LU factorisation of A
     for dof in bdofs
         A[dof, dof] = 1.0e60
     end
     flush!(A)
-    LUA = LUFactorization(A)
+    LUA = lu(A.cscmatrix)
 
     return MyPreconditionerLogPrimal{Tv, typeof(LUA)}(LUA, DA, bdofs, nmodes)
 end
 
-
-#@inline LinearAlgebra.ldiv!(x::AbstractArray, C::ExtendableSparseMatrix, b::AbstractArray) = x = C\b
 @inline LinearAlgebra.ldiv!(C::MyPreconditionerLogPrimal, b) = ldiv!(b, C, b)
 @inline function LinearAlgebra.ldiv!(y, C::MyPreconditionerLogPrimal{Tv, FAC}, b) where {Tv, FAC}
     #copyto!(y, b)
diff --git a/src/modelproblems/solvers_poisson_primal.jl b/src/modelproblems/solvers_poisson_primal.jl
index e73a200..543d496 100644
--- a/src/modelproblems/solvers_poisson_primal.jl
+++ b/src/modelproblems/solvers_poisson_primal.jl
@@ -38,13 +38,11 @@ function MyPreconditionerPrimal(A::ExtendableSparseMatrix{Tv, Ti}, bdofs, nmodes
         A[dof, dof] = 1.0e60
     end
     flush!(A)
-    LUA = LUFactorization(A)
+    LUA = lu(A.cscmatrix)
 
     return MyPreconditionerPrimal{Tv, typeof(LUA)}(LUA, DA, bdofs, nmodes)
 end
 
-
-#@inline LinearAlgebra.ldiv!(x::AbstractArray, C::ExtendableSparseMatrix, b::AbstractArray) = x = C\b
 @inline LinearAlgebra.ldiv!(C::MyPreconditionerPrimal, b) = ldiv!(b, C, b)
 @inline function LinearAlgebra.ldiv!(y, C::MyPreconditionerPrimal{Tv, FAC}, b) where {Tv, FAC}
     #copyto!(y, b)