diff --git a/sample/testNewton.cpp b/sample/testNewton.cpp
new file mode 100644
index 0000000..27415b8
--- /dev/null
+++ b/sample/testNewton.cpp
@@ -0,0 +1,27 @@
+#include <iostream>
+#include <cmath>
+
+using namespace std;
+
+#include "newton.hpp"
+
+using namespace CosmoTool;
+
+struct SimplePolynom
+{
+  double eval(double x) { return x*x*x+2*x-1; }
+  double derivative(double x) { return 3*x*x+2; }
+};
+
+int main()
+{
+  SimplePolynom poly;
+  double r;
+
+  double solution = -2*pow(2./(3*(9+sqrt(177.))), 1./3) + (pow(0.5*(9+sqrt(177.))/9, 1./3));
+
+  r = newtonSolver(3.0, poly);
+
+  cout << "Result = " << r << " delta = " << r-solution << endl;
+  return 0;
+}
diff --git a/src/newton.hpp b/src/newton.hpp
new file mode 100644
index 0000000..ca39804
--- /dev/null
+++ b/src/newton.hpp
@@ -0,0 +1,28 @@
+#ifndef _COSMOTOOL_NEWTON_HPP
+#define _COSMOTOOL_NEWTON_HPP
+
+namespace CosmoTool
+{
+  template<typename T, typename FunT>
+  T newtonSolver(T x0, FunT function, double residual = 1e-3)
+  {
+    T x, xold = x0;
+    T f_x = function.eval(x0);
+    T df_x = function.derivative(x0);
+    
+    x = xold - f_x/df_x;
+
+    while (abs(xold-x) > residual)
+      {	
+	xold = x;
+	f_x = function.eval(x);
+	df_x = function.derivative(x);
+	x = xold - f_x/df_x;
+      }
+	   
+    return x;
+  }
+
+};
+
+#endif