diff --git a/demo/mixed_poisson_dlrbnicsx/mesh_data/mesh.py b/demo/mixed_poisson_dlrbnicsx/mesh_data/mesh.py
index ef9f01c..cbcb46c 100644
--- a/demo/mixed_poisson_dlrbnicsx/mesh_data/mesh.py
+++ b/demo/mixed_poisson_dlrbnicsx/mesh_data/mesh.py
@@ -6,12 +6,12 @@
 gmsh.model.add("3d_mesh")
 
 gdim = 3
-lc = 5.e-2
+lc = 5.e-1
 
-gmsh.model.geo.addPoint(0., 0., 0., lc, 1)
-gmsh.model.geo.addPoint(1., 0., 0., lc, 2)
-gmsh.model.geo.addPoint(1., 1., 0., lc, 3)
-gmsh.model.geo.addPoint(0., 1., 0., lc, 4)
+gmsh.model.geo.addPoint(0., 0., 0., 0.2, 1)
+gmsh.model.geo.addPoint(1., 0., 0., 0.3, 2)
+gmsh.model.geo.addPoint(1., 1., 0., 0.4, 3)
+gmsh.model.geo.addPoint(0., 1., 0., 0.5, 4)
 
 gmsh.model.geo.addLine(1, 2, 1)
 gmsh.model.geo.addLine(2, 3, 2)
@@ -25,15 +25,15 @@
 
 # Create mesh
 # 8 = Frontal-Delaunay for Quads
-gmsh.option.setNumber("Mesh.Algorithm", 8)
+# gmsh.option.setNumber("Mesh.Algorithm", 8)
 # 2 = simple full-quad
-gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
+# gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
 # Apply recombination algorithm
-gmsh.option.setNumber("Mesh.RecombineAll", 1)
+# gmsh.option.setNumber("Mesh.RecombineAll", 1)
 # Mesh subdivision algorithm (1: all quadrangles)
-gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
+# gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
 # Minimum characteristic element size
-gmsh.option.setNumber("Mesh.MeshSizeMin", lc)
+gmsh.option.setNumber("Mesh.MeshSizeMin", 0.2)
 # Maximum characteristic element size
 gmsh.option.setNumber("Mesh.MeshSizeMax", lc)
 # Mesh generation
diff --git a/demo/time_dependent_diffusion/dolfinx_time_dependent_diffusion.py b/demo/time_dependent_diffusion/dolfinx_time_dependent_diffusion.py
index 7e98d48..01b4ef4 100644
--- a/demo/time_dependent_diffusion/dolfinx_time_dependent_diffusion.py
+++ b/demo/time_dependent_diffusion/dolfinx_time_dependent_diffusion.py
@@ -7,6 +7,12 @@
 from dolfinx.fem.petsc import \
     assemble_matrix, assemble_vector, create_vector, set_bc, apply_lifting
 
+import rbnicsx
+import rbnicsx.online
+import rbnicsx.backends
+
+import abc
+
 # Read mesh
 mesh_comm = MPI.COMM_WORLD
 gdim = 3
@@ -50,6 +56,9 @@
 xdmf.write_mesh(mesh)
 xdmf.write_function(uD_prev, 0)
 
+print("Set up snapshots matrix")
+snapshots_matrix = rbnicsx.backends.FunctionsList(VD)
+
 for i in range(num_steps):
     t += dt
     print(f"Time: {t}s")
@@ -61,6 +70,8 @@
     set_bc(b, bc)
     solver.solve(b, uD_sol.vector)
     uD_sol.x.scatter_forward()
+    
+    snapshots_matrix.append(uD_sol) # TODO see if uD_sol is overwritten every time
 
     uD_prev.x.array[:] = uD_sol.x.array
 
@@ -69,3 +80,72 @@
 
 xdmf.close()
 
+snapshots_matrix = rbnicsx.backends.FunctionsList(VD)
+
+class PODANNReducedProblem(abc.ABC):
+    def __init__(self, fem_space):
+        self._basis_functions = rbnicsx.backends.FunctionsList(fem_space)
+        u, v = ufl.TrialFunction(fem_space), ufl.TestFunction(fem_space)
+        self._inner_product = ufl.inner(u, v) * ufl.dx +\
+            ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
+        self._inner_product_action = \
+            rbnicsx.backends.bilinear_form_action(self._inner_product,
+                                                  part="real")
+
+# Maximum RB size
+Nmax = 30
+
+print(rbnicsx.io.TextBox("POD offline phase begins", fill="="))
+print("")
+
+print("Set up reduced problem")
+reduced_problem = PODANNReducedProblem(VD)
+
+print("")
+
+#for (mu_index, mu) in enumerate(training_set):
+    #print(rbnicsx.io.TextLine(str(mu_index+1), fill="#"))
+
+    #print("Parameter number ", (mu_index+1), "of", training_set.shape[0])
+    #print("High fidelity solve for mu =", mu)
+    #snapshot = problem_parametric.solve(mu)
+    #print(f"Solution array: {snapshot.x.array}")
+
+    #print("Update snapshots matrix")
+    #snapshots_matrix.append(snapshot)
+
+    #print("")
+
+print(rbnicsx.io.TextLine("Perform POD", fill="#"))
+eigenvalues, modes, _ = \
+    rbnicsx.backends.\
+    proper_orthogonal_decomposition(snapshots_matrix,
+                                    reduced_problem._inner_product_action,
+                                    N=Nmax, tol=1e-10)
+reduced_problem._basis_functions.extend(modes)
+reduced_size = len(reduced_problem._basis_functions)
+print("")
+
+print(rbnicsx.io.TextBox("POD-Galerkin offline phase ends", fill="="))
+
+positive_eigenvalues = np.where(eigenvalues > 0., eigenvalues, np.nan)
+singular_values = np.sqrt(positive_eigenvalues)
+
+print(f"Eigenvalues: {positive_eigenvalues}")
+
+plt.figure(figsize=[8, 10])
+xint = list()
+yval = list()
+
+for x, y in enumerate(eigenvalues[:len(reduced_problem._basis_functions)]):
+    yval.append(y)
+    xint.append(x+1)
+
+plt.plot(xint, yval, "*-", color="orange")
+plt.xlabel("Eigenvalue number", fontsize=18)
+plt.ylabel("Eigenvalue", fontsize=18)
+plt.xticks(xint)
+plt.yscale("log")
+plt.title("Eigenvalue decay", fontsize=24)
+plt.tight_layout()
+plt.savefig("eigenvalue_decay")