{% if (SSMethod == "kinsol") %}

  int fsolve (KINSysFn f, int N, N_Vector fval, void *user_data, std::string name)
  {
      int status;
      N_Vector u0, sc;
      void *kmem;

      u0 = N_VNew_Serial(N);
      N_VConst_Serial(0.0,u0); 
      N_VConst_Serial(0.0,fval); 

      sc = N_VNew_Serial(N);
      N_VConst_Serial(1.0,sc); 

      kmem = KINCreate();

      status = KINSetUserData (kmem, user_data);
      if (check_flag (&status, "KinSetUserData", 1)) throw KINSolverFailure (name, status);

      status = KINInit (kmem, f, u0);
      if (check_flag (&status, "KinInit", 1)) throw KINSolverFailure (name, status);

      status = KINDense (kmem, N);
      if (check_flag (&status, "KinDense", 1)) throw KINSolverFailure (name, status);

      status = KINSol (kmem, fval, KIN_NONE, sc, sc);
      if (check_flag (&status, "KINSol", 1)) throw KINSolverFailure (name, status);
      
      N_VDestroy_Serial(sc);
      N_VDestroy_Serial(u0);

      KINFree (&kmem);

      return 0;
  }

{% else %}

  int fsolve (int (*fss)(const gsl_vector *, void *user_data, gsl_vector *), int N, gsl_vector *fval, void *user_data, std::string name)
  {
      
      const gsl_multiroot_fsolver_type * T = gsl_multiroot_fsolver_hybrid;
      gsl_multiroot_fsolver * s = gsl_multiroot_fsolver_alloc (T, N);
      gsl_multiroot_function f = {fss, N, user_data};

      int status, iter;
      gsl_vector *x = gsl_vector_alloc (N);
      
      for (int i = 0; i < N; i++)
      {
         gsl_vector_set (x, i, 0.0);
      }

      gsl_multiroot_fsolver_set (s, &f, x);

      iter = 0;
      do
      {
         iter++;
         status = gsl_multiroot_fsolver_iterate (s);

         if ((status == GSL_EBADFUNC) ||
             (status == GSL_ENOPROG))
         {
            throw GSLSolverFailure(name, status);
         }
         
         status =  gsl_multiroot_test_residual (s->f, 1e-7);
      }
      while (status == GSL_CONTINUE && iter < 1000);

      for (int i = 0; i < N; i++)
      {
         gsl_vector_set (fval, i, gsl_vector_get (s->x, i));
      }

      gsl_vector_free (x);
      gsl_multiroot_fsolver_free (s);

      return 0;
  }

{% endif %}