diff --git a/examples/1_basic_workflow.jmd b/examples/1_basic_workflow.jmd
index 38bc883..5cadf9b 100644
--- a/examples/1_basic_workflow.jmd
+++ b/examples/1_basic_workflow.jmd
@@ -32,10 +32,13 @@ meta_p==meta_q
 ```
 
 ```julia
-g = ConScape.Grid(size(mov_prob)...,
-                  affinities=ConScape.graph_matrix_from_raster(mov_prob),
-                  qualities=hab_qual,
-                  costs=ConScape.MinusLog(), )
+g = ConScape.Grid(
+  size(mov_prob)...,
+  affinities=ConScape.graph_matrix_from_raster(mov_prob),
+  qualities=hab_qual,
+  costs=ConScape.MinusLog()
+)
+
 ```
 
 ```julia
@@ -58,16 +61,71 @@ g.target_qualities
 
 ```julia
 g.affinities
+g.affinities[1:100, 1:100]
+
 ```
 
 ```julia
 (g.nrows, g.ncols, g.nrows*g.ncols)
+g.costmatrix
 ```
 
 # Step 2: Habitat creation
 
 ```julia
-h = ConScape.GridRSP(g, θ=1.0);
+using ProfileView
+using BenchmarkTools
+
+_, targetnodes = idx_and_nodes = ConScape._targetidx_and_nodes(g)
+
+Z = Matrix{Float64}(undef, size(g.costmatrix, 1), length(idx_and_nodes[1]))
+
+B = Matrix(sparse(targetnodes,
+  1:length(targetnodes),
+  1.0,
+  size(g.costmatrix, 1),
+  length(targetnodes))
+)
+Z = LinearAlgebra.copy_similar(B, Float64)
+
+opt = ConScape.GridRSPoptimised(g; θ=1.0, idx_and_nodes, Z, B)
+orig = ConScape.GridRSPoriginal(g; θ=1.0)
+orig.Z .- opt.Z
+
+Pref = ConScape._Pref(g.affinities)
+W = ConScape._W(Pref, θ, g.costmatrix)
+A = SparseArrays.I - W
+F = lu(A)
+
+@allocated h = ConScape.GridRSPoptimised(g; θ=1.0, idx_and_nodes, Z, F)
+@allocated h = ConScape.GridRSPoptimised(g; θ=1.0)
+@allocated h = ConScape.GridRSPoriginal(g; θ=1.0)
+@b ConScape.GridRSPoptimised(g; θ=1.0, idx_and_nodes, Z, F)
+@b ConScape.GridRSPoptimised(g; θ=1.0)
+@b ConScape.GridRSPoriginal(g; θ=1.0)
+
+@profview ConScape.GridRSPoptimised(g; θ=1.0, idx_and_nodes, Z, F)
+@profview h = ConScape.GridRSPoriginal(g; θ=1.0)
+h = ConScape.GridRSPoptimised(g; θ=1.0)
+grsp = h
+h
+L = log.(h.Z)
+using GLMakie
+fig = Figure()
+a1 = Axis(fig[1, 1])
+a2 = Axis(fig[1, 2])
+A1 = broadcast(collect(A)) do x
+    isapprox(x, 0) ? NaN : x
+end
+Z1 = broadcast(h.Z) do x
+    isapprox(x, 0.0; atol=1e-5) ? NaN : x
+end
+p1 = Makie.heatmap!(a1, A1; colormap=:viridis)
+p2 = Makie.heatmap!(a2, Z1; colormap=:viridis)
+Makie.Colorbar(fig[1, 1, Right()], p1)
+Makie.Colorbar(fig[1, 2, Right()], p2)
+linkaxes!(a1, a2)
+
 ```
 
 show distance:
@@ -83,11 +141,11 @@ dists = ConScape.expected_cost(h);
 ```
 
 ```julia
-ConScape.plot_values(g, dists[:,4300], title="RSP Expected Cost Distance")
+ConScape.plot_values(g, dists[:, 4300], title="RSP Expected Cost Distance")
 ```
 
 ```julia
-ConScape.plot_values(g, dists[:,4300], title="RSP Expected Cost Distance")
+ConScape.plot_values(g, dists[:, 4300], title="RSP Expected Cost Distance")
 # savefig("figure_ECDistance.png")
 ```
 
diff --git a/src/gridrsp.jl b/src/gridrsp.jl
index 8e9f711..c5a2fb0 100644
--- a/src/gridrsp.jl
+++ b/src/gridrsp.jl
@@ -11,7 +11,59 @@ end
 
 Construct a GridRSP from a `g::Grid` based on the inverse temperature parameter `θ::Real`.
 """
-function GridRSP(g::Grid; θ=nothing)
+function GridRSPoptimised(g::Grid; 
+    θ=nothing,
+    idx_and_nodes=ConScape._targetidx_and_nodes(g),
+    targetidx=idx_and_nodes[1],
+    targetnodes=idx_and_nodes[2],
+    Z=Matrix(sparse(targetnodes,
+        1:length(targetnodes),
+        1.0,
+        size(g.costmatrix, 1),
+        length(targetnodes))),
+    F=nothing,
+)
+    Pref = ConScape._Pref(g.affinities)
+    W = ConScape._W(Pref, θ, g.costmatrix)
+    A = SparseArrays.I - W
+
+    @debug("Computing fundamental matrix of non-absorbing paths (Z). Please be patient...")
+
+    if isnothing(F)
+        F = lu(A)
+    else
+        lu!(F, A)
+    end
+    transposeoptype = SparseArrays.LibSuiteSparse.UMFPACK_A
+
+    # SparseArrays.UMFPACK._AqldivB_kernel!(Z, F, B, transposeoptype)
+
+    # This is basically SparseArrays.UMFPACK._AqldivB_kernel!
+    # But we unroll it to avoid copies or allocation of B
+    B = zeros(size(Z, 2))
+    i = 1
+    n = targetnodes[i]
+    for col in 1:size(Z, 2)
+        # Manually set the diagonal
+        if targetnodes[i] == col
+            B[col] = 1.0
+        end
+        SparseArrays.UMFPACK.solve!(view(Z, :, col), F, B, transposeoptype)
+        if targetnodes[i] == col
+            B[col] = 0.0
+            i += 1
+        end
+    end
+
+    # Check that values in Z are not too small:
+    if minimum(Z) * minimum(nonzeros(g.costmatrix .* W)) == 0
+        @warn "Warning: Z-matrix contains too small values, which can lead to inaccurate results! Check that the graph is connected or try decreasing θ."
+    end
+
+    return GridRSP(g, θ, Pref, W, Z)
+end
+
+function GridRSPoriginal(g::Grid; θ=nothing)
 
     Pref = _Pref(g.affinities)
     W    = _W(Pref, θ, g.costmatrix)