Supporting code from the Appendix

Kronecker Tensor Product of Matrices

SetAttributes[KroneckerTensor, OneIdentity]; 
    KroneckerTensor[u_?MatrixQ, v_?MatrixQ] := 
      Module[{w = Outer[Times, u, v]}, 
        Partition[
          Flatten[Transpose[
              w, {1, 3, 2, 4}]], Dimensions[
                w][[2]] Dimensions[
                w][[4]]]];
KroneckerTensor[u_, v_, w__] := 
        Fold[KroneckerTensor, u, {v, w}];
CircleTimes = KroneckerTensor;

Session p. 124 (alternatively to using MATLAB)

Matrices[
          n_] := Matrices[n] = (m = 2 n + 1; (Subscript[p, E]) = (1/4) + ϵ0; 
            (Subscript[p, W]) = (1/4) - ϵ0; (Subscript[p, N]) = (1/4); (Subscript[p, S]) = (1/4); 
            Id = SparseArray[{i_, i_} -> 1, {m, m}]; 
            PEW = SparseArray[{{i_, j_} /; j == i + 1 -> 
                    (Subscript[p, E]), {i_, j_} /; j == i - 1 -> (Subscript[p, W])}, {m, 
                  m}]; PNS = 
              SparseArray[{{i_, j_} /; j == i + 1 -> 
                    (Subscript[p, N]), {i_, j_} /; j == i - 1 -> (Subscript[p, S])}, {m, 
                  m}]; {PEW⊗Id + Id⊗PNS, 
              Id⊗Id});
\(ReturnProbability[ϵ_Real, 
        n_Integer] := \(({A, Id} = Matrices[n]; 
        Block[{ϵ0 = ϵ}, m = 2\ n + 1; ctr = n\ m + n + 1; 
          r = A[[All, ctr]]; 
          A[[All, ctr]] = 0; 
          q = LinearSolve[Id - A, r]; 
          q[[ctr]]])\); \({n == #, 
            N == \((2  # + 1)\)\^2 - 
                1, \(ϵ\_*\) == \((ϵ /. 
                  FindRoot[
                    ReturnProbability[ϵ, #] == 
                      1\/2, {ϵ, 0.06, 0.07}, 
                    WorkingPrecision -> MachinePrecision, 
                    PrecisionGoal -> 12\ ])\)} &\) /@ {10, 20, 40} // 
      TableForm\)