diff --git a/examples/algebraicSolver.m b/examples/algebraicSolver.m
index 15084f31b955b2fc67e5442a9a104c4c26881c7e..3ef04b1eb3a6f50650b74ff9ad1cf3380613f7b9 100644
--- a/examples/algebraicSolver.m
+++ b/examples/algebraicSolver.m
@@ -33,7 +33,7 @@ for p = 1:pmax
     % choose the iterative solver of choice: multigrid (with the variants
     % lowOrderVcycle and highOrderVcycle), pcg (with jacobi and additive
     % Schwarz preconditioner), and the cg solver
-    [solver, P] = chooseIterativeSolver(fes, blf, 'multigrid', 'lowOrderVcycle');
+    [solver, P] = IterativeSolver.choose(fes, blf, 'multigrid', 'lowOrderVcycle');
     solver.tol = 1e-6;
     solver.maxIter = 100;
 
diff --git a/examples/goafemIterativeSolver.m b/examples/goafemIterativeSolver.m
index 801bbfb0ba1caaa3660e4ee777db7969d4fcc842..50a7df3835d03592ed0f894e33081a62b8935177 100644
--- a/examples/goafemIterativeSolver.m
+++ b/examples/goafemIterativeSolver.m
@@ -37,7 +37,7 @@ for p = 1:pmax
     lfG.qrf = QuadratureRule.ofOrder(2*p);
 
     %% set up solver and lifting operator for nested iteration
-    [solver, P] = chooseIterativeSolver(fes, blf, 'multigrid', 'lowOrderVcycle');
+    [solver, P] = IterativeSolver.choose(fes, blf, 'multigrid', 'lowOrderVcycle');
     solver.tol = 1e-6;
     solver.maxIter = 1000;
 
diff --git a/lib/solvers/IterativeSolver.m b/lib/solvers/@IterativeSolver/IterativeSolver.m
similarity index 92%
rename from lib/solvers/IterativeSolver.m
rename to lib/solvers/@IterativeSolver/IterativeSolver.m
index 191676e17144f503abcdefd31db89cafbba62aec..35495b17807f3316424294e274af7cfb95f86356 100644
--- a/lib/solvers/IterativeSolver.m
+++ b/lib/solvers/@IterativeSolver/IterativeSolver.m
@@ -25,6 +25,10 @@
 %
 % tf = solver.isConverged(solver) returns state of convergence of the solver
 %   for each column of the right-hand side.
+%
+% [solver, P] = IterativeSolver.choose(fes, blf, class, variant) convenience
+%   factory method. Returns instances of IterativeSolver and Prolongation
+%   accoding to the given parameters.
 
 
 classdef IterativeSolver < handle
@@ -102,6 +106,11 @@ classdef IterativeSolver < handle
             obj.activeComponents = obj.activeComponents & ~tf;
         end
     end
+
+    %% convenience factory function
+    methods (Static, Access=public)
+        solver = choose(fes, blf, class, variant)
+    end
     
     %% validation functions to use within this class
     methods (Static, Access=protected)
diff --git a/lib/solvers/chooseIterativeSolver.m b/lib/solvers/@IterativeSolver/choose.m
similarity index 77%
rename from lib/solvers/chooseIterativeSolver.m
rename to lib/solvers/@IterativeSolver/choose.m
index 54bf5fbed9886c2f330ff2fae58b04ec04fba34f..ef5d4fa6dbc32592f5c496f354ff839080dbd3a7 100644
--- a/lib/solvers/chooseIterativeSolver.m
+++ b/lib/solvers/@IterativeSolver/choose.m
@@ -1,7 +1,18 @@
-% chooseIterativeSolver Return suitable solver (instance of subclass of
-%   IterativeSolver) and suitable Prolongation P.
+% choose Return suitable solver (instance of subclass of IterativeSolver) and
+%   suitable Prolongation P.
+%
+%   solver = choose(fes, blf, class) returns a solver of class 'class' for given
+%       finite element space fes and bilinear form blf. For each class, sensible
+%       default variants are chosen.
+%
+%   solver = choose(fes, blf, class, variant) returns a solver of class 'class'
+%       of specific variant 'variant'.
+%
+%   [solver, P] = choose(__) additionally to the solver, returns a suitable
+%       Prolongation P, that can be used to prolongate the solution of the
+%       solver to the next finer mesh.
 
-function [solver, P] = chooseIterativeSolver(fes, blf, class, variant)
+function [solver, P] = choose(fes, blf, class, variant)
 arguments
     fes FeSpace
     blf BilinearForm
@@ -48,7 +59,7 @@ switch class
         end
         solver = PcgSolver(preconditioner);
         
-    % geometric multigrid family
+    % geometric , Pmultigrid family
     case "multigrid"
         switch variant
             case {"", "lowOrderVcycle"}
@@ -78,5 +89,3 @@ switch class
 end
     
 end
-
-
diff --git a/runAdditiveSchwarzPreconditioner.m b/runAdditiveSchwarzPreconditioner.m
index 3e6d606bd1631a6eaad60923c6a2d12516e26ae9..971ffd8d1f8979c73c28c1216fafa76b00187343 100644
--- a/runAdditiveSchwarzPreconditioner.m
+++ b/runAdditiveSchwarzPreconditioner.m
@@ -32,10 +32,10 @@ function leveldata = runAdditiveSchwarzPreconditioner(pMax)
         lf.f = Constant(mesh, 1);
     
         % Choose iterative solver
-        solver = chooseIterativeSolver(fes, blf, "pcg", "additiveSchwarzLowOrder");
-        % solver = chooseIterativeSolver(fes, blf, "pcg", "additiveSchwarzHighOrder");
-        % solver = chooseIterativeSolver(fes, blf, "pcg", "iChol");
-        % solver = chooseIterativeSolver(fes, blf, "multigrid", "lowOrderVcycle");
+        solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzLowOrder");
+        % solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzHighOrder");
+        % solver = IterativeSolver.choose(fes, blf, "pcg", "iChol");
+        % solver = IterativeSolver.choose(fes, blf, "multigrid", "lowOrderVcycle");
     
         for k = 1:nLevels-1
             % Assemble matrices on intermediate meshes
diff --git a/runIterativeSolver.m b/runIterativeSolver.m
index 350c6b1f2801128319ed64ee03192f39facac566..ce2bcf8b6d8cd180066937d7c5843e07e05907f6 100644
--- a/runIterativeSolver.m
+++ b/runIterativeSolver.m
@@ -28,10 +28,10 @@ function leveldata = runIterativeSolver(maxNiter)
     lf.f = Constant(mesh, 1);
 
     % Choose iterative solver
-    solver = chooseIterativeSolver(fes, blf, "pcg", "additiveSchwarzLowOrder");
-    % solver = chooseIterativeSolver(fes, blf, "pcg", "additiveSchwarzHighOrder");
-    % solver = chooseIterativeSolver(fes, blf, "pcg", "iChol");
-    % solver = chooseIterativeSolver(fes, blf, "multigrid", "lowOrderVcycle");
+    solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzLowOrder");
+    % solver = IterativeSolver.choose(fes, blf, "pcg", "additiveSchwarzHighOrder");
+    % solver = IterativeSolver.choose(fes, blf, "pcg", "iChol");
+    % solver = IterativeSolver.choose(fes, blf, "multigrid", "lowOrderVcycle");
 
     % Initialize LevelData
     leveldata = LevelData('results');
diff --git a/tests/TestIterativeSolver.m b/tests/TestIterativeSolver.m
index cbde88457dc1aa0e1c8392b95ad3e36857882f06..6a32e6af2f72e24d59a80be2ceb5d94690187682 100644
--- a/tests/TestIterativeSolver.m
+++ b/tests/TestIterativeSolver.m
@@ -16,7 +16,7 @@ methods (Test)
     function firstStepDecreasesNorm(testCase, p, variant)
         [~, fes, blf, lf] = setupProblem(testCase, p);
         s = variant(1); v = variant(2);
-        solver = chooseIterativeSolver(fes, blf, s, v);
+        solver = IterativeSolver.choose(fes, blf, s, v);
         
         [xstar, A, F] = assembleData(testCase, blf, lf, fes);
         solver.setupSystemMatrix(A);
@@ -32,7 +32,7 @@ methods (Test)
     function laterStepDecreasesNorm(testCase, p, variant)
         [mesh, fes, blf, lf] = setupProblem(testCase, p);
         s = variant(1); v = variant(2);
-        solver = chooseIterativeSolver(fes, blf, s, v);
+        solver = IterativeSolver.choose(fes, blf, s, v);
         
         for k = 1:3
             [xstar, A, F] = assembleData(testCase, blf, lf, fes);
@@ -51,7 +51,7 @@ methods (Test)
     function solverIsLinear(testCase, p, variant)
         [mesh, fes, blf, lf] = setupProblem(testCase, p);
         s = variant(1); v = variant(2);
-        solver = chooseIterativeSolver(fes, blf, s, v);
+        solver = IterativeSolver.choose(fes, blf, s, v);
         
         for k = 1:3
             [xstar, A, F] = assembleData(testCase, blf, lf, fes);