diff --git a/MathLib/Nonlinear/Picard-impl.h b/MathLib/Nonlinear/Picard-impl.h
index 93b1174798a9559c33e93278d12ee1e13ed4760d..1052ba72b6ab422ed89528bab86b85b3ee719938 100644
--- a/MathLib/Nonlinear/Picard-impl.h
+++ b/MathLib/Nonlinear/Picard-impl.h
@@ -27,8 +27,8 @@ Picard::Picard()
 {
 }
 
-template<class F_PROBLEM, class T_VALUE>
-bool Picard::solve(F_PROBLEM &g,  const T_VALUE &x0, T_VALUE &x_new)
+template<class T_FUNCTOR, class T_VALUE>
+bool Picard::solve(T_FUNCTOR &functor,  const T_VALUE &x0, T_VALUE &x_new)
 {
     T_VALUE x_old(x0);
     T_VALUE dx(x0);
@@ -45,7 +45,7 @@ bool Picard::solve(F_PROBLEM &g,  const T_VALUE &x0, T_VALUE &x_new)
     double abs_error = -1.;
     double rel_error = -1.;
     for (itr_cnt=0; itr_cnt<_max_itr; itr_cnt++) {
-        g(x_old, x_new);
+        functor(x_old, x_new);
         dx = x_new;
         dx -= x_old;
 
diff --git a/MathLib/Nonlinear/Picard.h b/MathLib/Nonlinear/Picard.h
index 02f20e4e913c863e1fa419fd89f4bb13efb8e6e5..212e6e7efdfb58818d110113491f64285740ad15 100644
--- a/MathLib/Nonlinear/Picard.h
+++ b/MathLib/Nonlinear/Picard.h
@@ -49,13 +49,21 @@ public:
     /**
      * solve a nonlinear problem
      *
-     * \param g                 Fixed point function object (x=g(x))
+     * \tparam T_FUNCTOR        Function object type which supports an
+     * operator ()(\f$x_k\f$, \f$x_{k+1}\f$). \f$x_k\f$ is a previous iteration
+     * step value. \f$x_{k+1}\f$ is a new solution.
+     * \tparam T_VALUE          Data type of \f$x_k\f$ and \f$x_{k+1}\f$.
+     * Both scalar and vector types are available as far as the following conditions
+     * are satisfied
+     * - T_VALUE has a copy constructor (for non-primitive data types)
+     * - MathLib::norm(T_VALUE) exists
+     * \param functor           Fixed point function object (\f$x_{k+1}=g(x_k)\f$)
      * \param x0                Initial guess
      * \param x_new             Solution
      * \return true if converged
      */
-    template<class F_PROBLEM, class T_VALUE>
-    bool solve(F_PROBLEM &g,  const T_VALUE &x0, T_VALUE &x_new);
+    template<class T_FUNCTOR, class T_VALUE>
+    bool solve(T_FUNCTOR &functor,  const T_VALUE &x0, T_VALUE &x_new);
 
     /// return the number of iterations
     std::size_t getNIterations() const {return _n_iterations; }