diff --git a/examples/5_anisotropic_square/anisotropic_square.py b/examples/5_anisotropic_square/anisotropic_square.py
new file mode 100644
index 0000000..db86b94
--- /dev/null
+++ b/examples/5_anisotropic_square/anisotropic_square.py
@@ -0,0 +1,113 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from itertools import product
+from anisotropic_square_engine import (AnisotropicSquareEngine as engine,
+                                       AnisotropicSquareEnginePoles as plsEx)
+from rrompy.sampling.standard import SamplingEngineStandard as SES
+from rrompy.reduction_methods.standard import NearestNeighbor as NN
+from rrompy.reduction_methods.pivoted.greedy import (
+                               RationalInterpolantGreedyPivotedGreedy as RIGPG)
+from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
+from rrompy.parameter import localSparseGrid as LSG
+
+zs, Ls = [10., 50.], [.2, 1.2]
+z0, L0, n = np.mean(zs), np.mean(Ls), 50
+murange = [[zs[0], Ls[0]], [zs[-1], Ls[-1]]]
+np.random.seed(4020)
+mu = [zs[0] + np.random.rand() * (zs[-1] - zs[0]),
+      Ls[0] + np.random.rand() * (Ls[-1] - Ls[0])]
+solver = engine(z0, L0, n)
+
+fighandles = []
+params = {"POD": True, "nTestPoints": 100, "greedyTol": 1e-4, "S": 3,
+          "polybasisMarginal": "MONOMIAL_GAUSSIAN",
+          "polybasis": "LEGENDRE", 'samplerPivot':QS(zs, "UNIFORM"),
+          'trainSetGenerator':QS(zs, "UNIFORM"),
+          'errorEstimatorKind':"INTERPOLATORY", "MMarginal": 2, "SMarginal": 3,
+          "greedyTolMarginal": 2.5e-2, "samplerMarginalGrid":LSG(Ls),
+          "radialDirectionalWeightsMarginal": [2.], "matchingWeight": 1.}
+
+for tol, kind in product([1., 3.], ["HARD", "SOFT"]):
+    print("Testing cutoff tolerance {} with kind {}.".format(tol, kind))
+    params['cutOffTolerance'] = tol
+    params['cutOffKind'] = kind
+
+    approx = RIGPG([0], solver, mu0 = [z0, L0], approx_state = True,
+                   approxParameters = params, verbosity = 5)
+    approx.setupApprox("LAST")
+    
+    print("--- Approximant ---")
+    approx.plotApprox(mu, name = 'u_app')
+    approx.plotHF(mu, name = 'u_HF')
+    approx.plotErr(mu, name = 'err_app')
+    approx.plotRes(mu, name = 'res_app')
+    normErr = approx.normErr(mu)[0]
+    normSol = approx.normHF(mu)[0]
+    normRes = approx.normRes(mu)[0]
+    normRHS = approx.normRHS(mu)[0]
+    print("SolNorm:\t{:.5e}\nErr_app: \t{:.5e}\nErrRel_app:\t{:.5e}".format(
+                                          normSol, normErr, normErr / normSol))
+    print("RHSNorm:\t{:.5e}\nRes_app: \t{:.5e}\nResRel_app:\t{:.5e}".format(
+                                          normRHS, normRes, normRes / normRHS))
+    
+    print("--- Closest snapshot ---")
+    eng = SES(solver, verbosity = 0)
+    eng.nsamples = approx.samplingEngine.nsamplesCoalesced
+    eng.mus = approx.samplingEngine.musCoalesced
+    eng.samples = approx.samplingEngine.samples_fullCoalesced
+    paramsNN = {'S':eng.nsamples, 'POD':False,
+                'sampler':QS(murange, "UNIFORM")}
+    approxNN = NN(solver, mu0 = [z0, L0], approx_state = True,
+                  approxParameters = paramsNN, verbosity = 0)
+    approxNN.setSamples(eng)
+    approxNN.plotApprox(mu, name = 'u_close')
+    approxNN.plotHF(mu, name = 'u_HF')
+    approxNN.plotErr(mu, name = 'err_close')
+    approxNN.plotRes(mu, name = 'res_close')
+    normErr = approxNN.normErr(mu)[0]
+    normSol = approxNN.normHF(mu)[0]
+    normRes = approxNN.normRes(mu)[0]
+    normRHS = approxNN.normRHS(mu)[0]
+    print("SolNorm:\t{:.5e}\nErr_close:\t{:.5e}\nErrRel_close:\t{:.5e}".format(
+                                          normSol, normErr, normErr / normSol))
+    print("RHSNorm:\t{:.5e}\nRes_close:\t{:.5e}\nResRel_close:\t{:.5e}".format(
+                                          normRHS, normRes, normRes / normRHS))
+    
+    verb = approx.verbosity
+    approx.verbosity = 0
+    tspace = np.linspace(Ls[0], Ls[-1], 100)
+    for j, t in enumerate(tspace):
+        plsE = plsEx(t, 0., zs[-1])
+        pls = approx.getPoles([None, t])
+        pls[np.abs(np.imag(pls)) > 1e-5] = np.nan
+        if j == 0:
+            polesE = np.empty((len(tspace), len(plsE)))
+            poles = np.empty((len(tspace), len(pls)))
+            polesE[:] = np.nan
+        if len(plsE) > polesE.shape[1]:
+            nanR = np.empty((len(tspace), len(plsE) - polesE.shape[1]))
+            nanR[:] = np.nan
+            polesE = np.hstack((polesE, nanR))
+        polesE[j, : len(plsE)] = np.real(plsE)
+        poles[j] = np.real(pls)
+    approx.verbosity = verb
+    fighandles += [plt.figure(figsize = (17, 5))]
+    ax1 = fighandles[-1].add_subplot(1, 2, 1)
+    ax2 = fighandles[-1].add_subplot(1, 2, 2)
+    ax1.plot(poles, tspace)
+    ax1.set_ylim(Ls)
+    ax1.set_xlabel('mu_1')
+    ax1.set_ylabel('mu_2')
+    ax1.grid()
+    ax2.plot(polesE, tspace, 'k-.', linewidth = 1)
+    ax2.plot(poles, tspace)
+    for mm in approx.musMarginal:
+        ax2.plot(zs, [mm[0, 0]] * 2, 'k--', linewidth = 1)
+    ax2.set_xlim(zs)
+    ax2.set_ylim(Ls)
+    ax2.set_xlabel('mu_1')
+    ax2.set_ylabel('mu_2')
+    ax2.grid()
+    plt.show()
+
+    print("\n")
diff --git a/examples/5_anisotropic_square/anisotropic_square_engine.py b/examples/5_anisotropic_square/anisotropic_square_engine.py
new file mode 100644
index 0000000..56eb191
--- /dev/null
+++ b/examples/5_anisotropic_square/anisotropic_square_engine.py
@@ -0,0 +1,67 @@
+import numpy as np
+import fenics as fen
+import mshr
+import ufl
+from rrompy.hfengines.linear_problem import HelmholtzProblemEngine
+from rrompy.solver.fenics import fenONE, fenZERO, fenics2Sparse
+
+class AnisotropicSquareEngine(HelmholtzProblemEngine):
+    def __init__(self, k2:float, L2:float, n:int):
+        super().__init__(mu0 = [k2, L2])
+        self._affinePoly = True
+        self.nAs = 3
+        self.npar = 2
+        self.rescalingExp = [1., 1.]
+        mesh = fen.UnitSquareMesh(n, n)
+        mesh = mshr.generate_mesh(mshr.Rectangle(fen.Point(0., 0.),
+                                                 fen.Point(1., 1.)), n)
+        self.V = fen.FunctionSpace(mesh, "P", 1)
+        eps = 1e-6
+        self.DirichletBoundary = lambda x, on_boundary: (on_boundary
+                                                         and x[1] < eps)
+        self.NeumannBoundary = "REST"
+        x, y = fen.SpatialCoordinate(mesh)[:]
+        self.NeumannDatum = ufl.conditional(ufl.ge(y, 1. - eps),
+                                            fen.cos(np.pi * x), fenZERO)
+        self.forcingTerm = ufl.conditional(ufl.ge(y, .5), fenONE, fenZERO) * (
+                  5 * ufl.conditional(ufl.lt(x, .1), fenONE, fenZERO)
+                - 5 * ufl.conditional(ufl.And(ufl.gt(x, .2), ufl.lt(x, .3)),
+                                      fenONE, fenZERO)
+                + 10 * ufl.conditional(ufl.And(ufl.gt(x, .45), ufl.lt(x, .55)),
+                                       fenONE, fenZERO)
+                - 5 * ufl.conditional(ufl.And(ufl.gt(x, .7), ufl.lt(x, .8)),
+                                      fenONE, fenZERO)
+                + 5 * ufl.conditional(ufl.gt(x, .9), fenONE, fenZERO))
+    
+    def buildA(self):
+        """Build terms of operator of linear system."""
+        if self.thAs[0] is None: self.thAs = self.getMonomialWeights(self.nAs)
+        if self.As[0] is None:
+            self.autoSetDS()
+            DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
+                                           self.DirichletBoundary)
+            a0 = fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx
+            self.As[0] = fenics2Sparse(a0, {}, DirichletBC0, 1)
+        if self.As[1] is None:
+            self.autoSetDS()
+            DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
+                                           self.DirichletBoundary)
+            a1 = - fen.dot(self.u, self.v) * fen.dx
+            self.As[1] = fenics2Sparse(a1, {}, DirichletBC0, 0)
+        if self.As[2] is None:
+            self.autoSetDS()
+            DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
+                                           self.DirichletBoundary)
+            a2 = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
+            self.As[2] = fenics2Sparse(a2, {}, DirichletBC0, 0)
+
+def AnisotropicSquareEnginePoles(L2:float, k2min:float, k2max:float):
+    poles = []
+    for alpha in np.arange(np.ceil((k2max) ** .5 / np.pi)):
+        p = (np.pi * alpha) ** 2.
+        pkmin = np.ceil(max(0., (k2min - p) * 4 / L2) ** .5 / np.pi)
+        pkmin += 1 - (pkmin % 2)
+        pkmax = np.floor(max(0., (k2max - p) * 4 / L2) ** .5 / np.pi)
+        for beta in np.arange(pkmin, pkmax + 1, 2):
+            poles += [p + L2 * (np.pi * beta / 2.) ** 2.]
+    return np.unique(poles)