diff --git a/festim/__init__.py b/festim/__init__.py
index bd29c6078..e17c850c1 100644
--- a/festim/__init__.py
+++ b/festim/__init__.py
@@ -102,6 +102,7 @@
 from .concentration.concentration import Concentration
 from .initial_condition import InitialCondition
 from .concentration.mobile import Mobile
+from .nonlinear_problem import Problem
 from .concentration.theta import Theta
 
 from .concentration.traps.trap import Trap
diff --git a/festim/concentration/traps/extrinsic_trap.py b/festim/concentration/traps/extrinsic_trap.py
index 3b446a621..ebda566fa 100644
--- a/festim/concentration/traps/extrinsic_trap.py
+++ b/festim/concentration/traps/extrinsic_trap.py
@@ -1,4 +1,6 @@
 from festim import Trap, as_constant_or_expression
+from fenics import NewtonSolver, MPI
+import warnings
 
 
 class ExtrinsicTrapBase(Trap):
@@ -14,6 +16,7 @@ def __init__(
         relative_tolerance=1e-10,
         maximum_iterations=30,
         linear_solver=None,
+        preconditioner=None,
         **kwargs,
     ):
         """Inits ExtrinsicTrap
@@ -36,18 +39,48 @@ def __init__(
                 If None, the default fenics linear solver will be used ("umfpack").
                 More information can be found at: https://fenicsproject.org/pub/tutorial/html/._ftut1017.html.
                 Defaults to None.
+            preconditioner (str, optional): preconditioning method for the newton solver,
+                options can be veiwed by print(list_krylov_solver_preconditioners()).
+                Defaults to "default".
         """
         super().__init__(k_0, E_k, p_0, E_p, materials, density=None, id=id)
         self.absolute_tolerance = absolute_tolerance
         self.relative_tolerance = relative_tolerance
         self.maximum_iterations = maximum_iterations
         self.linear_solver = linear_solver
+        self.preconditioner = preconditioner
 
+        self.newton_solver = None
         for name, val in kwargs.items():
             setattr(self, name, as_constant_or_expression(val))
         self.density_previous_solution = None
         self.density_test_function = None
 
+    @property
+    def newton_solver(self):
+        return self._newton_solver
+
+    @newton_solver.setter
+    def newton_solver(self, value):
+        if value is None:
+            self._newton_solver = value
+        elif isinstance(value, NewtonSolver):
+            if self._newton_solver:
+                print("Settings for the Newton solver will be overwritten")
+            self._newton_solver = value
+        else:
+            raise TypeError("accepted type for newton_solver is fenics.NewtonSolver")
+
+    def define_newton_solver(self):
+        """Creates the Newton solver and sets its parameters"""
+        self.newton_solver = NewtonSolver(MPI.comm_world)
+        self.newton_solver.parameters["error_on_nonconvergence"] = False
+        self.newton_solver.parameters["absolute_tolerance"] = self.absolute_tolerance
+        self.newton_solver.parameters["relative_tolerance"] = self.relative_tolerance
+        self.newton_solver.parameters["maximum_iterations"] = self.maximum_iterations
+        self.newton_solver.parameters["linear_solver"] = self.linear_solver
+        self.newton_solver.parameters["preconditioner"] = self.preconditioner
+
 
 class ExtrinsicTrap(ExtrinsicTrapBase):
     """
diff --git a/festim/concentration/traps/traps.py b/festim/concentration/traps/traps.py
index f91a5880d..a19ee18b8 100644
--- a/festim/concentration/traps/traps.py
+++ b/festim/concentration/traps/traps.py
@@ -109,26 +109,23 @@ def define_variational_problem_extrinsic_traps(self, dx, dt, T):
                 self.extrinsic_formulations.append(trap.form_density)
         self.sub_expressions.extend(expressions_extrinsic)
 
+    def define_newton_solver_extrinsic_traps(self):
+        for trap in self:
+            if isinstance(trap, festim.ExtrinsicTrapBase):
+                trap.define_newton_solver()
+
     def solve_extrinsic_traps(self):
         for trap in self:
             if isinstance(trap, festim.ExtrinsicTrapBase):
                 du_t = f.TrialFunction(trap.density[0].function_space())
                 J_t = f.derivative(trap.form_density, trap.density[0], du_t)
-                problem = f.NonlinearVariationalProblem(
-                    trap.form_density, trap.density[0], [], J_t
-                )
-                solver = f.NonlinearVariationalSolver(problem)
-                solver.parameters["newton_solver"][
-                    "absolute_tolerance"
-                ] = trap.absolute_tolerance
-                solver.parameters["newton_solver"][
-                    "relative_tolerance"
-                ] = trap.relative_tolerance
-                solver.parameters["newton_solver"][
-                    "maximum_iterations"
-                ] = trap.maximum_iterations
-                solver.parameters["newton_solver"]["linear_solver"] = trap.linear_solver
-                solver.solve()
+                problem = festim.Problem(J_t, trap.form_density, [])
+
+                f.begin(
+                    "Solving nonlinear variational problem."
+                )  # Add message to fenics logs
+                trap.newton_solver.solve(problem, trap.density[0].vector())
+                f.end()
 
     def update_extrinsic_traps_density(self):
         for trap in self:
diff --git a/festim/h_transport_problem.py b/festim/h_transport_problem.py
index eea9c7722..f973e3b75 100644
--- a/festim/h_transport_problem.py
+++ b/festim/h_transport_problem.py
@@ -1,5 +1,6 @@
 from fenics import *
 import festim
+import warnings
 
 
 class HTransportProblem:
@@ -22,6 +23,7 @@ class HTransportProblem:
             ct2, ...)
         v (fenics.TestFunction): the test function
         u_n (fenics.Function): the "previous" function
+        newton_solver (fenics.NewtonSolver): Newton solver for solving the nonlinear problem
         bcs (list): list of fenics.DirichletBC for H transport
     """
 
@@ -36,6 +38,7 @@ def __init__(self, mobile, traps, T, settings, initial_conditions) -> None:
         self.u = None
         self.v = None
         self.u_n = None
+        self.newton_solver = None
 
         self.boundary_conditions = []
         self.bcs = None
@@ -43,6 +46,21 @@ def __init__(self, mobile, traps, T, settings, initial_conditions) -> None:
         self.V_CG1 = None
         self.expressions = []
 
+    @property
+    def newton_solver(self):
+        return self._newton_solver
+
+    @newton_solver.setter
+    def newton_solver(self, value):
+        if value is None:
+            self._newton_solver = value
+        elif isinstance(value, NewtonSolver):
+            if self._newton_solver:
+                print("Settings for the Newton solver will be overwritten")
+            self._newton_solver = value
+        else:
+            raise TypeError("accepted type for newton_solver is fenics.NewtonSolver")
+
     def initialise(self, mesh, materials, dt=None):
         """Assigns BCs, create suitable function space, initialise
         concentration fields, define variational problem
@@ -59,6 +77,7 @@ def initialise(self, mesh, materials, dt=None):
             self.mobile.volume_markers = mesh.volume_markers
             self.mobile.T = self.T
         self.attribute_flux_boundary_conditions()
+
         # Define functions
         self.define_function_space(mesh)
         self.initialise_concentrations()
@@ -71,12 +90,14 @@ def initialise(self, mesh, materials, dt=None):
             self.mobile.create_form_post_processing(self.V_DG1, materials, mesh.dx)
 
         self.define_variational_problem(materials, mesh, dt)
+        self.define_newton_solver()
 
         # Boundary conditions
         print("Defining boundary conditions")
         self.create_dirichlet_bcs(materials, mesh)
         if self.settings.transient:
             self.traps.define_variational_problem_extrinsic_traps(mesh.dx, dt, self.T)
+            self.traps.define_newton_solver_extrinsic_traps()
 
     def define_function_space(self, mesh):
         """Creates a suitable function space for H transport problem
@@ -201,6 +222,22 @@ def define_variational_problem(self, materials, mesh, dt=None):
         self.F = F
         self.expressions = expressions
 
+    def define_newton_solver(self):
+        """Creates the Newton solver and sets its parameters"""
+        self.newton_solver = NewtonSolver(MPI.comm_world)
+        self.newton_solver.parameters["error_on_nonconvergence"] = False
+        self.newton_solver.parameters["absolute_tolerance"] = (
+            self.settings.absolute_tolerance
+        )
+        self.newton_solver.parameters["relative_tolerance"] = (
+            self.settings.relative_tolerance
+        )
+        self.newton_solver.parameters["maximum_iterations"] = (
+            self.settings.maximum_iterations
+        )
+        self.newton_solver.parameters["linear_solver"] = self.settings.linear_solver
+        self.newton_solver.parameters["preconditioner"] = self.settings.preconditioner
+
     def attribute_flux_boundary_conditions(self):
         """Iterates through self.boundary_conditions, checks if it's a FluxBC
         and its field is 0, and assign fluxes to self.mobile
@@ -239,7 +276,6 @@ def update(self, t, dt):
             t (float): the current time (s)
             dt (festim.Stepsize): the stepsize
         """
-
         festim.update_expressions(self.expressions, t)
 
         converged = False
@@ -264,28 +300,16 @@ def solve_once(self):
             int, bool: number of iterations for reaching convergence, True if
                 converged else False
         """
-
         if self.J is None:  # Define the Jacobian
             du = TrialFunction(self.u.function_space())
             J = derivative(self.F, self.u, du)
         else:
             J = self.J
-        problem = NonlinearVariationalProblem(self.F, self.u, self.bcs, J)
-        solver = NonlinearVariationalSolver(problem)
-        solver.parameters["newton_solver"]["error_on_nonconvergence"] = False
-        solver.parameters["newton_solver"][
-            "absolute_tolerance"
-        ] = self.settings.absolute_tolerance
-        solver.parameters["newton_solver"][
-            "relative_tolerance"
-        ] = self.settings.relative_tolerance
-        solver.parameters["newton_solver"][
-            "maximum_iterations"
-        ] = self.settings.maximum_iterations
-        solver.parameters["newton_solver"][
-            "linear_solver"
-        ] = self.settings.linear_solver
-        nb_it, converged = solver.solve()
+        problem = festim.Problem(J, self.F, self.bcs)
+
+        begin("Solving nonlinear variational problem.")  # Add message to fenics logs
+        nb_it, converged = self.newton_solver.solve(problem, self.u.vector())
+        end()
 
         return nb_it, converged
 
diff --git a/festim/nonlinear_problem.py b/festim/nonlinear_problem.py
new file mode 100644
index 000000000..f06bbd47e
--- /dev/null
+++ b/festim/nonlinear_problem.py
@@ -0,0 +1,31 @@
+import fenics as f
+
+
+class Problem(f.NonlinearProblem):
+    """
+    Class to set up a nonlinear variational problem (F(u, v)=0) to solve
+    by the Newton method based on the form of the variational problem, the Jacobian
+    form of the variational problem, and the boundary conditions
+
+    Args:
+        J (ufl.Form): the Jacobian form of the variational problem
+        F (ufl.Form): the form of the variational problem
+        bcs (list): list of fenics.DirichletBC
+    """
+
+    def __init__(self, J, F, bcs):
+        self.jacobian_form = J
+        self.residual_form = F
+        self.bcs = bcs
+        self.assembler = f.SystemAssembler(
+            self.jacobian_form, self.residual_form, self.bcs
+        )
+        f.NonlinearProblem.__init__(self)
+
+    def F(self, b, x):
+        """Assembles the RHS in Ax=b and applies the boundary conditions"""
+        self.assembler.assemble(b, x)
+
+    def J(self, A, x):
+        """Assembles the LHS in Ax=b and applies the boundary conditions"""
+        self.assembler.assemble(A)
diff --git a/festim/settings.py b/festim/settings.py
index 182ca7a49..533fba778 100644
--- a/festim/settings.py
+++ b/festim/settings.py
@@ -22,9 +22,12 @@ class Settings:
             the formulation will be computed only once at the beggining.
             Else it will be computed at each time step. Defaults to True.
         linear_solver (str, optional): linear solver method for the newton solver,
-            options can be veiwed by print(list_linear_solver_methods()).
+            options can be viewed by print(list_linear_solver_methods()).
             More information can be found at: https://fenicsproject.org/pub/tutorial/html/._ftut1017.html.
             Defaults to None, for the newton solver this is: "umfpack".
+        preconditioner (str, optional): preconditioning method for the newton solver,
+            options can be viewed by print(list_krylov_solver_preconditioners()).
+            Defaults to "default".
 
     Attributes:
         transient (bool): transient or steady state sim
@@ -40,6 +43,7 @@ class Settings:
         traps_element_type (str): Finite element used for traps.
         update_jacobian (bool):
         linear_solver (str): linear solver method for the newton solver
+        precondtitioner (str): preconditioning method for the newton solver
     """
 
     def __init__(
@@ -54,6 +58,7 @@ def __init__(
         traps_element_type="CG",
         update_jacobian=True,
         linear_solver=None,
+        preconditioner="default",
     ):
         # TODO maybe transient and final_time are redundant
         self.transient = transient
@@ -68,3 +73,4 @@ def __init__(
         self.traps_element_type = traps_element_type
         self.update_jacobian = update_jacobian
         self.linear_solver = linear_solver
+        self.preconditioner = preconditioner
diff --git a/festim/temperature/temperature_solver.py b/festim/temperature/temperature_solver.py
index bb1d09d66..54eea2130 100644
--- a/festim/temperature/temperature_solver.py
+++ b/festim/temperature/temperature_solver.py
@@ -1,6 +1,7 @@
 import festim
 import fenics as f
 import sympy as sp
+import warnings
 
 
 class HeatTransferProblem(festim.Temperature):
@@ -21,10 +22,14 @@ class HeatTransferProblem(festim.Temperature):
             If None, the default fenics linear solver will be used ("umfpack").
             More information can be found at: https://fenicsproject.org/pub/tutorial/html/._ftut1017.html.
             Defaults to None.
+        preconditioner (str, optional): preconditioning method for the newton solver,
+            options can be veiwed by print(list_krylov_solver_preconditioners()).
+            Defaults to "default".
 
     Attributes:
         F (fenics.Form): the variational form of the heat transfer problem
         v_T (fenics.TestFunction): the test function
+        newton_solver (fenics.NewtonSolver): Newton solver for solving the nonlinear problem
         initial_condition (festim.InitialCondition): the initial condition
         sub_expressions (list): contains time dependent fenics.Expression to
             be updated
@@ -41,6 +46,7 @@ def __init__(
         relative_tolerance=1e-10,
         maximum_iterations=30,
         linear_solver=None,
+        preconditioner="default",
     ) -> None:
         super().__init__()
         self.transient = transient
@@ -49,12 +55,29 @@ def __init__(
         self.relative_tolerance = relative_tolerance
         self.maximum_iterations = maximum_iterations
         self.linear_solver = linear_solver
+        self.preconditioner = preconditioner
 
         self.F = 0
         self.v_T = None
         self.sources = []
         self.boundary_conditions = []
         self.sub_expressions = []
+        self.newton_solver = None
+
+    @property
+    def newton_solver(self):
+        return self._newton_solver
+
+    @newton_solver.setter
+    def newton_solver(self, value):
+        if value is None:
+            self._newton_solver = value
+        elif isinstance(value, f.NewtonSolver):
+            if self._newton_solver:
+                print("Settings for the Newton solver will be overwritten")
+            self._newton_solver = value
+        else:
+            raise TypeError("accepted type for newton_solver is fenics.NewtonSolver")
 
     @property
     def initial_condition(self):
@@ -105,20 +128,21 @@ def create_functions(self, materials, mesh, dt=None):
         self.define_variational_problem(materials, mesh, dt)
         self.create_dirichlet_bcs(mesh.surface_markers)
 
+        if not self.newton_solver:
+            self.define_newton_solver()
+
         if not self.transient:
             print("Solving stationary heat equation")
             dT = f.TrialFunction(self.T.function_space())
             JT = f.derivative(self.F, self.T, dT)
-            problem = f.NonlinearVariationalProblem(
-                self.F, self.T, self.dirichlet_bcs, JT
-            )
-            solver = f.NonlinearVariationalSolver(problem)
-            newton_solver_prm = solver.parameters["newton_solver"]
-            newton_solver_prm["absolute_tolerance"] = self.absolute_tolerance
-            newton_solver_prm["relative_tolerance"] = self.relative_tolerance
-            newton_solver_prm["maximum_iterations"] = self.maximum_iterations
-            newton_solver_prm["linear_solver"] = self.linear_solver
-            solver.solve()
+            problem = festim.Problem(JT, self.F, self.dirichlet_bcs)
+
+            f.begin(
+                "Solving nonlinear variational problem."
+            )  # Add message to fenics logs
+            self.newton_solver.solve(problem, self.T.vector())
+            f.end()
+
             self.T_n.assign(self.T)
 
     def define_variational_problem(self, materials, mesh, dt=None):
@@ -197,6 +221,16 @@ def define_variational_problem(self, materials, mesh, dt=None):
                 for surf in bc.surfaces:
                     self.F += -bc.form * self.v_T * mesh.ds(surf)
 
+    def define_newton_solver(self):
+        """Creates the Newton solver and sets its parameters"""
+        self.newton_solver = f.NewtonSolver(f.MPI.comm_world)
+        self.newton_solver.parameters["error_on_nonconvergence"] = False
+        self.newton_solver.parameters["absolute_tolerance"] = self.absolute_tolerance
+        self.newton_solver.parameters["relative_tolerance"] = self.relative_tolerance
+        self.newton_solver.parameters["maximum_iterations"] = self.maximum_iterations
+        self.newton_solver.parameters["linear_solver"] = self.linear_solver
+        self.newton_solver.parameters["preconditioner"] = self.preconditioner
+
     def create_dirichlet_bcs(self, surface_markers):
         """Creates a list of fenics.DirichletBC and add time dependent
         expressions to .sub_expressions
@@ -228,16 +262,14 @@ def update(self, t):
             # Solve heat transfers
             dT = f.TrialFunction(self.T.function_space())
             JT = f.derivative(self.F, self.T, dT)  # Define the Jacobian
-            problem = f.NonlinearVariationalProblem(
-                self.F, self.T, self.dirichlet_bcs, JT
-            )
-            solver = f.NonlinearVariationalSolver(problem)
-            newton_solver_prm = solver.parameters["newton_solver"]
-            newton_solver_prm["absolute_tolerance"] = self.absolute_tolerance
-            newton_solver_prm["relative_tolerance"] = self.relative_tolerance
-            newton_solver_prm["maximum_iterations"] = self.maximum_iterations
-            newton_solver_prm["linear_solver"] = self.linear_solver
-            solver.solve()
+            problem = festim.Problem(JT, self.F, self.dirichlet_bcs)
+
+            f.begin(
+                "Solving nonlinear variational problem."
+            )  # Add message to fenics logs
+            self.newton_solver.solve(problem, self.T.vector())
+            f.end()
+
             self.T_n.assign(self.T)
 
     def is_steady_state(self):
diff --git a/test/h_transport_problem/test_solving.py b/test/h_transport_problem/test_solving.py
index a292b7374..ee4a3442d 100644
--- a/test/h_transport_problem/test_solving.py
+++ b/test/h_transport_problem/test_solving.py
@@ -1,5 +1,6 @@
 import festim
 import fenics as f
+import pytest
 
 
 def test_default_dt_min_value():
@@ -28,6 +29,7 @@ def test_default_dt_min_value():
     my_problem = festim.HTransportProblem(
         festim.Mobile(), festim.Traps([]), festim.Temperature(200), my_settings, []
     )
+    my_problem.define_newton_solver()
     my_problem.u = f.Function(V)
     my_problem.u_n = f.Function(V)
     my_problem.v = f.TestFunction(V)
@@ -52,6 +54,7 @@ def test_solve_once_jacobian_is_none():
     my_problem = festim.HTransportProblem(
         festim.Mobile(), festim.Traps([]), festim.Temperature(200), my_settings, []
     )
+    my_problem.define_newton_solver()
     my_problem.u = f.Function(V)
     my_problem.u_n = f.Function(V)
     my_problem.v = f.TestFunction(V)
@@ -79,6 +82,7 @@ def test_solve_once_returns_false():
     my_problem = festim.HTransportProblem(
         festim.Mobile(), festim.Traps([]), festim.Temperature(200), my_settings, []
     )
+    my_problem.define_newton_solver()
     my_problem.u = f.Function(V)
     my_problem.u_n = f.Function(V)
     my_problem.v = f.TestFunction(V)
@@ -94,8 +98,15 @@ def test_solve_once_returns_false():
     assert not converged
 
 
-def test_solve_once_linear_solver_mumps():
-    """Checks that solve_once() works when an alternative linear solver is used rather than the default"""
+@pytest.mark.parametrize("preconditioner", ["default", "icc"])
+def test_solve_once_linear_solver_gmres(preconditioner):
+    """
+    Checks that solve_once() works when an alternative linear solver is used
+    with/without a preconditioner rather than the default
+
+    Args:
+        preconditioner (str): the preconditioning method
+    """
     # build
     mesh = f.UnitIntervalMesh(8)
     V = f.FunctionSpace(mesh, "CG", 1)
@@ -104,11 +115,13 @@ def test_solve_once_linear_solver_mumps():
         absolute_tolerance=1e-10,
         relative_tolerance=1e-10,
         maximum_iterations=50,
-        linear_solver="mumps",
+        linear_solver="gmres",
+        preconditioner=preconditioner,
     )
     my_problem = festim.HTransportProblem(
         festim.Mobile(), festim.Traps([]), festim.Temperature(200), my_settings, []
     )
+    my_problem.define_newton_solver()
     my_problem.u = f.Function(V)
     my_problem.u_n = f.Function(V)
     my_problem.v = f.TestFunction(V)
@@ -122,3 +135,56 @@ def test_solve_once_linear_solver_mumps():
 
     # test
     assert converged
+
+
+class Test_solve_once_with_custom_solver:
+    """
+    Checks that a custom newton sovler can be used
+    """
+
+    def sim(self):
+        """Defines a model"""
+        mesh = f.UnitIntervalMesh(8)
+        V = f.FunctionSpace(mesh, "CG", 1)
+
+        my_settings = festim.Settings(
+            absolute_tolerance=1e-10,
+            relative_tolerance=1e-10,
+            maximum_iterations=50,
+        )
+        my_problem = festim.HTransportProblem(
+            festim.Mobile(), festim.Traps([]), festim.Temperature(200), my_settings, []
+        )
+        my_problem.define_newton_solver()
+        my_problem.u = f.Function(V)
+        my_problem.u_n = f.Function(V)
+        my_problem.v = f.TestFunction(V)
+        my_problem.F = (
+            (my_problem.u - my_problem.u_n) * my_problem.v * f.dx
+            + 1 * my_problem.v * f.dx
+            + f.dot(f.grad(my_problem.u), f.grad(my_problem.v)) * f.dx
+        )
+        return my_problem
+
+    def test_custom_solver(self):
+        """Solves the system using the built-in solver and using the f.NewtonSolver"""
+
+        # solve with the built-in solver
+        problem_1 = self.sim()
+        problem_1.solve_once()
+
+        # solve with the custom solver
+        problem_2 = self.sim()
+        problem_2.newton_solver = f.NewtonSolver()
+        problem_2.newton_solver.parameters["absolute_tolerance"] = (
+            problem_1.settings.absolute_tolerance
+        )
+        problem_2.newton_solver.parameters["relative_tolerance"] = (
+            problem_1.settings.relative_tolerance
+        )
+        problem_2.newton_solver.parameters["maximum_iterations"] = (
+            problem_1.settings.maximum_iterations
+        )
+        problem_2.solve_once()
+
+        assert (problem_1.u.vector() == problem_2.u.vector()).all()
diff --git a/test/heat_transfer_problem/test_solving.py b/test/heat_transfer_problem/test_solving.py
index f5cca7179..73fe404d9 100644
--- a/test/heat_transfer_problem/test_solving.py
+++ b/test/heat_transfer_problem/test_solving.py
@@ -3,9 +3,15 @@
 import fenics as f
 
 
-def test_create_functions_linear_solver_mumps():
-    """Checks that the function created by create_functions() has the expected value when an
-    alternative linear solver is used rather than the default"""
+@pytest.mark.parametrize("preconditioner", ["default", "icc"])
+def test_create_functions_linear_solver_gmres(preconditioner):
+    """
+    Checks that the function created by create_functions() has the expected value when an
+    alternative linear solver is used with/without a preconditioner rather than the default
+
+    Args:
+        preconditioner (str): the preconditioning method
+    """
 
     mesh = festim.MeshFromRefinements(10, size=0.1)
 
@@ -21,7 +27,8 @@ def test_create_functions_linear_solver_mumps():
         absolute_tolerance=1e-03,
         relative_tolerance=1e-10,
         maximum_iterations=30,
-        linear_solver="mumps",
+        linear_solver="gmres",
+        preconditioner=preconditioner,
     )
     my_problem.boundary_conditions = bcs
 
@@ -29,3 +36,52 @@ def test_create_functions_linear_solver_mumps():
     my_problem.create_functions(materials=materials, mesh=mesh)
 
     assert my_problem.T(0.05) == pytest.approx(1)
+
+
+class Test_solve_once_with_custom_solver:
+    """
+    Checks that a custom newton sovler can be used
+    """
+
+    def sim(self):
+        """Defines a model"""
+        bcs = [
+            festim.DirichletBC(surfaces=[1, 2], value=1, field="T"),
+        ]
+
+        my_problem = festim.HeatTransferProblem(
+            transient=False,
+            absolute_tolerance=1e-03,
+            relative_tolerance=1e-10,
+            maximum_iterations=30,
+        )
+        my_problem.boundary_conditions = bcs
+        return my_problem
+
+    def test_custom_solver(self):
+        """Solves the system using the built-in solver and using the f.NewtonSolver"""
+        mesh = festim.MeshFromRefinements(10, size=0.1)
+        materials = festim.Materials(
+            [festim.Material(id=1, D_0=1, E_D=0, thermal_cond=1)]
+        )
+        mesh.define_measures(materials)
+
+        # solve with the built-in solver
+        problem_1 = self.sim()
+        problem_1.create_functions(materials=materials, mesh=mesh)
+
+        # solve with the custom solver
+        problem_2 = self.sim()
+        problem_2.newton_solver = f.NewtonSolver()
+        problem_2.newton_solver.parameters["absolute_tolerance"] = (
+            problem_1.absolute_tolerance
+        )
+        problem_2.newton_solver.parameters["relative_tolerance"] = (
+            problem_1.relative_tolerance
+        )
+        problem_2.newton_solver.parameters["maximum_iterations"] = (
+            problem_1.maximum_iterations
+        )
+        problem_2.create_functions(materials=materials, mesh=mesh)
+
+        assert (problem_1.T.vector() == problem_2.T.vector()).all()
diff --git a/test/system/test_misc.py b/test/system/test_misc.py
index 0c0728248..2b4cecaec 100644
--- a/test/system/test_misc.py
+++ b/test/system/test_misc.py
@@ -1,4 +1,5 @@
 import festim as F
+import fenics as f
 import numpy as np
 import pytest
 import os
@@ -333,6 +334,173 @@ def test_materials_setter():
     assert my_model.materials is test_materials
 
 
+class TestFestimProblem:
+    """Tests the methods of the festim.Problem class"""
+
+    # define the FESTIM problem
+    mesh = f.UnitIntervalMesh(8)
+    V = f.FunctionSpace(mesh, "CG", 1)
+    u = f.Function(V)
+    v = f.TestFunction(V)
+    s = u * v * f.dx + v * f.dx
+    J = f.derivative(s, u)
+    problem = F.Problem(J, s, [])
+
+    # define the fenics assembler used in festim.Problem
+    assembler = f.SystemAssembler(J, s, [])
+    x = f.PETScVector()
+
+    def test_F(self):
+        """
+        Creates two epty matrices and checks
+        that the festim.Problem.J properly assembles the RHS of Ax=b
+        """
+        b1 = f.PETScVector()
+        b2 = f.PETScVector()
+
+        self.problem.F(b1, self.x)
+        self.assembler.assemble(b2, self.x)
+
+        assert (b1 == b2).all()
+
+    def test_J(self):
+        """
+        Creates two epty matrices and checks
+        that the festim.Problem.J properly assembles the LHS of Ax=b
+        """
+        A1 = f.PETScMatrix()
+        A2 = f.PETScMatrix()
+
+        self.problem.J(A1, self.x)
+        self.assembler.assemble(A2)
+
+        assert (A1.array() == A2.array()).all()
+
+
+class TestWarningsCustomSolver:
+    """
+    Creates a simulation object and checks that a TypeError (UserWarning) is raised
+    when the newton_solver attribute is given a value of the wrong type (solver is overwritten)
+    """
+
+    def sim(self):
+        """Defines a model"""
+        my_sim = F.Simulation()
+        my_sim.mesh = F.MeshFromVertices([1, 2, 3])
+        my_mat = F.Materials(
+            [F.Material(id=1, D_0=1, E_D=0, thermal_cond=1, heat_capacity=1, rho=1)]
+        )
+        my_sim.materials = my_mat
+        my_sim.T = F.HeatTransferProblem(
+            transient=True, initial_condition=F.InitialCondition(field="T", value=1)
+        )
+        my_sim.traps = F.ExtrinsicTrap(
+            1,
+            1,
+            1,
+            1,
+            my_mat,
+            phi_0=1,
+            n_amax=2,
+            n_bmax=2,
+            eta_a=3,
+            eta_b=4,
+            f_a=5,
+            f_b=6,
+        )
+        # add source to the HeatTransferProblem
+        my_sim.settings = F.Settings(
+            transient=True,
+            absolute_tolerance=1e-10,
+            relative_tolerance=1e-10,
+            final_time=1,
+        )
+        my_sim.dt = F.Stepsize(1)
+        return my_sim
+
+    @pytest.mark.parametrize("value", ["coucou", [0, 0], -1.0])
+    def test_wrong_type_solver_h_transport(self, value):
+        """
+        Checks that a TypeError is raised when the newton_solver attribute
+        of the HTransportProblem class is given a value of the wrong type
+        """
+        problem = self.sim()
+        problem.initialise()
+
+        with pytest.raises(
+            TypeError,
+            match="accepted type for newton_solver is fenics.NewtonSolver",
+        ):
+            problem.h_transport_problem.newton_solver = value
+
+    @pytest.mark.parametrize("value", ["coucou", [0, 0], -1.0])
+    def test_wrong_type_solver_heat_transport(self, value):
+        """
+        Checks that a TypeError is raised when the newton_solver attribute
+        of the HeatTransferProblem class is given a value of the wrong type
+        """
+        problem = self.sim()
+
+        with pytest.raises(
+            TypeError,
+            match="accepted type for newton_solver is fenics.NewtonSolver",
+        ):
+            problem.T.newton_solver = value
+
+    @pytest.mark.parametrize("value", ["coucou", [0, 0], -1.0])
+    def test_wrong_type_solver_ex_trap(self, value):
+        """
+        Checks that a TypeError is raised when the newton_solver attribute
+        of the ExtrinsicTrap class is given a value of the wrong type
+        """
+        problem = self.sim()
+
+        with pytest.raises(
+            TypeError,
+            match="accepted type for newton_solver is fenics.NewtonSolver",
+        ):
+            problem.traps[0].newton_solver = value
+
+    def test_warn_solver_h_transport(self, capsys):
+        """
+        Checks that a warning is raised when the newton_solver attribute
+        of the HTransportProblem class is overwritten
+        """
+        problem = self.sim()
+        problem.initialise()
+        problem.h_transport_problem.newton_solver = f.NewtonSolver()
+        assert (
+            "Settings for the Newton solver will be overwritten"
+            in capsys.readouterr().out
+        )
+
+    def test_warn_solver_heat_transport(self, capsys):
+        """
+        Checks that a warning is raised when the newton_solver attribute
+        of the HeatTransferProblem class is overwritten
+        """
+        problem = self.sim()
+        problem.initialise()
+        problem.T.newton_solver = f.NewtonSolver()
+        assert (
+            "Settings for the Newton solver will be overwritten"
+            in capsys.readouterr().out
+        )
+
+    def test_warn_solver_ex_trap(self, capsys):
+        """
+        Checks that a warning is raised when the newton_solver attribute
+        of the ExtrinsicTrap class is overwritten
+        """
+        problem = self.sim()
+        problem.initialise()
+        problem.traps[0].newton_solver = f.NewtonSolver()
+        assert (
+            "Settings for the Newton solver will be overwritten"
+            in capsys.readouterr().out
+        )
+
+
 def test_error_raised_when_no_IC_heat_transfer():
     """
     Checks that an error is raised when no initial condition is provided for
diff --git a/test/unit/test_traps/test_extrinsic_trap.py b/test/unit/test_traps/test_extrinsic_trap.py
index e36ba4010..5eb44c92f 100644
--- a/test/unit/test_traps/test_extrinsic_trap.py
+++ b/test/unit/test_traps/test_extrinsic_trap.py
@@ -87,3 +87,5 @@ def test_solver_parameters(self):
         self.my_trap.relative_tolerance = 1
         self.my_trap.maximum_iterations = 1
         self.my_trap.linear_solver = "mumps"
+        self.my_trap.preconditioner = "icc"
+        self.my_trap.newton_solver = f.NewtonSolver()
diff --git a/test/unit/test_traps/test_neutron_induced_trap.py b/test/unit/test_traps/test_neutron_induced_trap.py
index 302d81a6f..1d6a441ac 100644
--- a/test/unit/test_traps/test_neutron_induced_trap.py
+++ b/test/unit/test_traps/test_neutron_induced_trap.py
@@ -146,7 +146,8 @@ class TestNeutronInducedTrapSolverParameters:
         absolute_tolerance=2.5,
         relative_tolerance=1.2e-10,
         maximum_iterations=13,
-        linear_solver="mumps",
+        linear_solver="gmres",
+        preconditioner="icc",
     )
 
     def test_attributes_from_instanciation(self):
@@ -158,7 +159,8 @@ def test_attributes_from_instanciation(self):
         assert self.my_trap.absolute_tolerance == 2.5
         assert self.my_trap.relative_tolerance == 1.2e-10
         assert self.my_trap.maximum_iterations == 13
-        assert self.my_trap.linear_solver == "mumps"
+        assert self.my_trap.linear_solver == "gmres"
+        assert self.my_trap.preconditioner == "icc"
 
     def test_attributes_change_since_instanciation(self):
         """