diff --git a/CHANGELOG.md b/CHANGELOG.md
index 23855d71..2be087a5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,9 +8,15 @@
 
   - new example on coupled Stokes-Darcy (Example264)
 
+### Added
+
+  - example `Example301` now also demonstrates a nonlinear problem solved by an iterative linear solver
+    with preconditioning.
+
 ### Changed
 
   - `solve` uses now the residual equation for the linear systems
+  - facelift `Example250`.
 
 ### Fixed
 
diff --git a/Project.toml b/Project.toml
index b4a0aa63..60666180 100644
--- a/Project.toml
+++ b/Project.toml
@@ -33,6 +33,7 @@ ExtendableGrids = "1.10.3"
 ExtendableSparse = "1.5.3"
 ForwardDiff = "0.10.35,1"
 GridVisualize = "1.8.1"
+IncompleteLU = "0.2.1"
 LinearAlgebra = "1.9"
 LinearSolve = "2, 3"
 Metis = "1.5.0"
@@ -55,6 +56,7 @@ julia = "1.9"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
+IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
 OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
@@ -65,4 +67,4 @@ TetGen = "c5d3f3f7-f850-59f6-8a2e-ffc6dc1317ea"
 Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
 
 [targets]
-test = ["Aqua", "ExampleJuggler", "ExplicitImports", "Metis", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "SimplexGridFactory", "StaticArrays", "Test", "TetGen", "Triangulate"]
+test = ["Aqua", "ExampleJuggler", "ExplicitImports", "IncompleteLU", "Metis", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "SimplexGridFactory", "StaticArrays", "Test", "TetGen", "Triangulate"]
diff --git a/docs/Project.toml b/docs/Project.toml
index 10e6594a..e36b6994 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -8,6 +8,7 @@ ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
 ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
+IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
diff --git a/examples/Example250_NSELidDrivenCavity.jl b/examples/Example250_NSELidDrivenCavity.jl
index 840937fc..1491235a 100644
--- a/examples/Example250_NSELidDrivenCavity.jl
+++ b/examples/Example250_NSELidDrivenCavity.jl
@@ -21,44 +21,74 @@ The computed solution for the default parameters looks like this:
 module Example250_NSELidDrivenCavity
 
 using ExtendableFEM
-using GridVisualize
+using Triangulate
 using ExtendableGrids
+using SimplexGridFactory
 using LinearAlgebra
+using GridVisualize
 using Test #hide
 
-function kernel_nonlinear!(result, u_ops, qpinfo)
-    u, ∇u, p = view(u_ops, 1:2), view(u_ops, 3:6), view(u_ops, 7)
+
+function create_cone(h)
+    builder = SimplexGridBuilder(; Generator = Triangulate)
+
+    ## points
+    p1 = point!(builder, -1, 0)
+    p2 = point!(builder, 1, 0)
+    p3 = point!(builder, 0, -2)
+
+    ## top face
+    facetregion!(builder, 1)
+    facet!(builder, p1, p2)
+
+    ## other faces
+    facetregion!(builder, 2)
+    facet!(builder, p2, p3)
+    facet!(builder, p3, p1)
+
+    cellregion!(builder, 1)
+    maxvolume!(builder, h)
+    regionpoint!(builder, 0, -0.5)
+
+    return simplexgrid(builder)
+end
+
+const 𝕀 = [1 0; 0 1]
+
+function NSE_kernel!(result, u_ops, qpinfo)
+
+    u = tensor_view(u_ops, 1, TDVector(2))
+    v = tensor_view(result, 1, TDVector(2))
+    ∇u = tensor_view(u_ops, 3, TDMatrix(2))
+    ∇v = tensor_view(result, 3, TDMatrix(2))
+    p = tensor_view(u_ops, 7, TDScalar())
+    q = tensor_view(result, 7, TDScalar())
     μ = qpinfo.params[1]
-    result[1] = dot(u, view(∇u, 1:2))
-    result[2] = dot(u, view(∇u, 3:4))
-    result[3] = μ * ∇u[1] - p[1]
-    result[4] = μ * ∇u[2]
-    result[5] = μ * ∇u[3]
-    result[6] = μ * ∇u[4] - p[1]
-    result[7] = -(∇u[1] + ∇u[4])
+
+    tmul!(v, ∇u, u)
+    ∇v .= μ .* ∇u .- p[1] .* 𝕀
+    q[1] = -dot(∇u, 𝕀)
+
     return nothing
 end
 
-function boundarydata!(result, qpinfo)
+
+## boundary function
+function u_boundary!(result, qpinfo)
     result[1] = 1
     result[2] = 0
     return nothing
 end
 
-function initialgrid_cone()
-    xgrid = ExtendableGrid{Float64, Int32}()
-    xgrid[Coordinates] = Array{Float64, 2}([-1 0; 0 -2; 1 0]')
-    xgrid[CellNodes] = Array{Int32, 2}([1 2 3]')
-    xgrid[CellGeometries] = VectorOfConstants{ElementGeometries, Int}(Triangle2D, 1)
-    xgrid[CellRegions] = ones(Int32, 1)
-    xgrid[BFaceRegions] = Array{Int32, 1}([1, 2, 3])
-    xgrid[BFaceNodes] = Array{Int32, 2}([1 2; 2 3; 3 1]')
-    xgrid[BFaceGeometries] = VectorOfConstants{ElementGeometries, Int}(Edge1D, 3)
-    xgrid[CoordinateSystem] = Cartesian2D
-    return xgrid
-end
 
-function main(; μ_final = 0.0005, order = 2, nrefs = 5, Plotter = nothing, kwargs...)
+function main(;
+        μ_final = 0.0005,        # flow parameter
+        order = 2,               # FE order of the flow field (pressure order is order-1)
+        h = 1.0e-3,              # grid cell volume
+        nrefs = 1,               # additional grid refinements
+        Plotter = nothing,
+        kwargs...
+    )
 
     ## prepare parameter field
     extra_params = Array{Float64, 1}([max(μ_final, 0.05)])
@@ -70,17 +100,20 @@ function main(; μ_final = 0.0005, order = 2, nrefs = 5, Plotter = nothing, kwar
 
     assign_unknown!(PD, u)
     assign_unknown!(PD, p)
-    assign_operator!(PD, NonlinearOperator(kernel_nonlinear!, [id(u), grad(u), id(p)]; params = extra_params, kwargs...))
-    assign_operator!(PD, InterpolateBoundaryData(u, boundarydata!; regions = 3))
-    assign_operator!(PD, HomogeneousBoundaryData(u; regions = [1, 2]))
-    assign_operator!(PD, FixDofs(p; dofs = [1]))
+    assign_operator!(PD, NonlinearOperator(NSE_kernel!, [id(u), grad(u), id(p)]; params = extra_params, kwargs...))
 
+    ## boundary data
+    assign_operator!(PD, InterpolateBoundaryData(u, u_boundary!; regions = 1))
+    assign_operator!(PD, HomogeneousBoundaryData(u; regions = [2]))
 
     ## grid
-    xgrid = uniform_refine(initialgrid_cone(), nrefs)
+    xgrid = uniform_refine(create_cone(h), nrefs)
 
     ## prepare FESpace
-    FES = [FESpace{H1Pk{2, 2, order}}(xgrid), FESpace{H1Pk{1, 2, order - 1}}(xgrid)]
+    FES = [
+        FESpace{H1Pk{2, 2, order}}(xgrid),
+        FESpace{H1Pk{1, 2, order - 1}}(xgrid),
+    ]
 
     ## prepare plots
     plt = GridVisualizer(; Plotter = Plotter, layout = (1, 2), clear = true, size = (1600, 800))
@@ -93,7 +126,16 @@ function main(; μ_final = 0.0005, order = 2, nrefs = 5, Plotter = nothing, kwar
     while (true)
         step += 1
         @info "Step $step : solving for μ=$(extra_params[1])"
-        sol, SC = ExtendableFEM.solve(PD, FES, SC; return_config = true, target_residual = 1.0e-10, maxiterations = 20, kwargs...)
+        sol, SC = ExtendableFEM.solve(
+            PD,
+            FES,
+            SC;
+            return_config = true,
+            target_residual = 1.0e-10,
+            maxiterations = 20,
+            kwargs...
+        )
+
         if step == 1
             initialize!(PE, sol)
         end
@@ -116,9 +158,10 @@ function main(; μ_final = 0.0005, order = 2, nrefs = 5, Plotter = nothing, kwar
 end
 
 generateplots = ExtendableFEM.default_generateplots(Example250_NSELidDrivenCavity, "example250.png") #hide
-function runtests() #hide
-    sol, plt = main(; nrefs = 3, μ_final = 0.005) #hide
-    @test sum(view(sol[1])) ≈ 9.501630403050289 #hide
-    return nothing #hide
-end #hide
+function runtests()                                                                                  #hide
+    sol, plt = main(; μ_final = 0.005)                                                               #hide
+    sum(view(sol[1])) ≈ 237.24628017878518                                                           #hide
+    return nothing                                                                                   #hide
+end                                                                                                  #hide
+
 end # module
diff --git a/examples/Example301_PoissonProblem.jl b/examples/Example301_PoissonProblem.jl
index 17609303..78cdde64 100644
--- a/examples/Example301_PoissonProblem.jl
+++ b/examples/Example301_PoissonProblem.jl
@@ -15,46 +15,99 @@ parameters looks like this:
 
 ![](example301.png)
 
+This examples uses an iterative solver with an IncompleteLU preconditioner as the default solver.
+It can be changed via the arguments 'method_linear' and 'precon_linear', see the runtests function
+for more examples.
+
 =#
 
 module Example301_PoissonProblem
 
 using ExtendableFEM
 using ExtendableGrids
-using Test #hide
+using LinearSolve
+using IncompleteLU
+using LinearAlgebra
+using Test
+
+function f!(result, qpinfo)
+    result[1] = qpinfo.params[1] * (1.7^2 + 3.9^2) * sin(1.7 * qpinfo.x[1]) * cos(3.9 * qpinfo.x[2])
+    return nothing
+end
 
-function f!(fval, qpinfo)
-    fval[1] = qpinfo.x[1] * qpinfo.x[2] * qpinfo.x[3]
+function u!(result, qpinfo)
+    result[1] = sin(1.7 * qpinfo.x[1]) * cos(3.9 * qpinfo.x[2])
     return nothing
 end
 
-function main(; μ = 1.0, nrefs = 3, Plotter = nothing, kwargs...)
+function main(;
+        μ = 1.0,
+        nrefs = 4,
+        method_linear = KrylovJL_GMRES(),
+        precon_linear = method_linear == KrylovJL_GMRES() ? IncompleteLU.ilu : nothing,
+        Plotter = nothing,
+        kwargs...
+    )
 
     ## problem description
     PD = ProblemDescription()
     u = Unknown("u"; name = "potential")
     assign_unknown!(PD, u)
-    assign_operator!(PD, BilinearOperator([grad(u)]; factor = μ, kwargs...))
-    assign_operator!(PD, LinearOperator(f!, [id(u)]; kwargs...))
-    assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4))
+    assign_operator!(PD, BilinearOperator([grad(u)]; factor = μ))
+    assign_operator!(PD, LinearOperator(f!, [id(u)]; params = [μ]))
+    assign_operator!(PD, InterpolateBoundaryData(u, u!; regions = 1:6))
 
     ## discretize
     xgrid = uniform_refine(grid_unitcube(Tetrahedron3D), nrefs)
     FES = FESpace{H1P2{1, 3}}(xgrid)
 
     ## solve
-    sol = solve(PD, FES; kwargs...)
+    sol = solve(PD, FES; method_linear, precon_linear, kwargs...)
+
+    ## compute error
+    function exact_error!(result, u, qpinfo)
+        u!(result, qpinfo)
+        result .-= u
+        result .= result .^ 2
+        return nothing
+    end
+    ErrorIntegratorExact = ItemIntegrator(exact_error!, [id(u)]; quadorder = 8)
+
+    ## calculate error
+    error = evaluate(ErrorIntegratorExact, sol)
+    L2error = sqrt(sum(view(error, 1, :)))
+    @info "L2 error = $L2error"
 
     ## plot
     plt = plot([id(u)], sol; Plotter = Plotter)
 
-    return sol, plt
+    return L2error, plt
 end
 
-generateplots = ExtendableFEM.default_generateplots(Example301_PoissonProblem, "example301.png") #hide
-function runtests() #hide
-    sol, plt = main() #hide
-    @test sum(sol.entries) ≈ 21.874305144549524 #hide
-    return nothing #hide
-end #hide
+generateplots = ExtendableFEM.default_generateplots(Example301_PoissonProblem, "example301.png")
+function runtests()
+    expected_error = 8.56e-5
+
+    ## test direct solver
+    L2error, plt = main(; nrefs = 4)
+    @test L2error <= expected_error
+
+    ## test iterative solver with IncompleteLU (fastest)
+    method_linear = KrylovJL_GMRES()
+    precon_linear = IncompleteLU.ilu
+    L2error, plt = main(; method_linear, precon_linear, nrefs = 4)
+    @test L2error <= expected_error
+
+    ## test other iterative solvers
+    method_linear = KrylovJL_GMRES(precs = (A, p) -> (Diagonal(A), I))
+    precon_linear = nothing
+    L2error, plt = main(; method_linear, precon_linear, nrefs = 4)
+    @test L2error <= expected_error
+
+    ## also working:
+    ## method_linear = KrylovJL_GMRES(precs = (A, p) -> (AMGCLWrap.AMGPrecon(ruge_stuben(A)), I))
+    ## method_linear = KrylovJL_GMRES(precs = (A, p) -> (AlgebraicMultigrid.aspreconditioner(ruge_stuben(A)), I))
+
+    return nothing
+end
 end # module