Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 24054)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 24055)
@@ -162,4 +162,5 @@
 					./shared/Numerics/cubic.cpp\
 					./shared/Numerics/NewtonSolveDnorm.cpp\
+					./shared/Numerics/ODE1.cpp\
 					./shared/Numerics/extrema.cpp\
 					./shared/Numerics/legendre.cpp\
Index: /issm/trunk-jpl/src/c/shared/Numerics/ODE1.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/ODE1.cpp	(revision 24055)
+++ /issm/trunk-jpl/src/c/shared/Numerics/ODE1.cpp	(revision 24055)
@@ -0,0 +1,23 @@
+#include <math.h>
+#include "./types.h"
+#include "../Exceptions/exceptions.h"
+
+IssmDouble ODE1(IssmDouble alpha,IssmDouble beta,IssmDouble Si, IssmDouble dt,int method){
+	/* solve the following equation:
+	 *
+	 *   dS/dt = alpha S + beta
+	 *
+	 *   method 0: Forward Euler (explicit)
+	 *   method 1: backward Euler (implicit)
+	 *   method 2: Crank Nicolson
+	 *
+	 *   return S^{i+1} based on  Si, dt, alpha and beta
+	 */
+
+	switch(method){
+		case 0: return Si*(1.+alpha*dt) + beta*dt;
+		case 1: return (Si+beta*dt)/(1.-alpha*dt);
+		case 2: return (Si*(1.+alpha*dt/2.) + beta*dt)/(1-alpha*dt/2.)
+		default: _error_("not supported yet");
+	}
+}
Index: /issm/trunk-jpl/src/c/shared/Numerics/numerics.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/numerics.h	(revision 24054)
+++ /issm/trunk-jpl/src/c/shared/Numerics/numerics.h	(revision 24055)
@@ -41,4 +41,5 @@
 
 int         NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n,IssmDouble dnorm);
+IssmDouble  ODE1(IssmDouble alpha,IssmDouble beta,IssmDouble Si, IssmDouble dt,int method);
 
 #endif //ifndef _NUMERICS_H_