4 files changed, 98 insertions, 5 deletions

diff --git a/packages/sundials/src/Test.hs b/packages/sundials/src/Main.hs
index a99582a..9e8bc63 100644
--- a/packages/sundials/src/Test.hs
+++ b/packages/sundials/src/Main.hs
@@ -150,15 +150,45 @@ C.include "helpers.h"
 -- }
 
 main = do
-  res <- [C.block| int { sunindextype NEQ = 1;        /* number of dependent vars. */
-                         N_Vector y = NULL;           /* empty vector for storing solution */
-                         void *arkode_mem = NULL;     /* empty ARKode memory structure */
+  res <- [C.block| int { /* general problem variables */
+                         int flag;                       /* reusable error-checking flag */
+                         N_Vector y = NULL;              /* empty vector for storing solution */
+                         SUNMatrix A = NULL;             /* empty matrix for linear solver */
+                         SUNLinearSolver LS = NULL;      /* empty linear solver object */
+                         void *arkode_mem = NULL;        /* empty ARKode memory structure */
+                         FILE *UFID;
+                         realtype t, tout;
+                         long int nst, nst_a, nfe, nfi, nsetups, nje, nfeLS, nni, ncfn, netf;
+
+                         /* general problem parameters */
+                         realtype T0 = RCONST(0.0);    /* initial time */
+                         realtype Tf = RCONST(10.0);   /* final time */
+                         realtype dTout = RCONST(1.0); /* time between outputs */
+                         sunindextype NEQ = 1;         /* number of dependent vars. */
+                         realtype reltol = 1.0e-6;     /* tolerances */
+                         realtype abstol = 1.0e-10;
+                         realtype lamda  = -100.0;     /* stiffness parameter */
+  
+                         /* Initial diagnostics output */
+                         printf("\nAnalytical ODE test problem:\n");
+                         printf("    lamda = %"GSYM"\n",    lamda);
+                         printf("   reltol = %.1"ESYM"\n",  reltol);
+                         printf("   abstol = %.1"ESYM"\n\n",abstol);
+
+                         /* Initialize data structures */
                          y = N_VNew_Serial(NEQ);      /* Create serial vector for solution */
                          if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
-
                          N_VConst(0.0, y);            /* Specify initial condition */
                          arkode_mem = ARKodeCreate(); /* Create the solver memory */
                          if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
+
+                         /* Call ARKodeInit to initialize the integrator memory and specify the */
+                         /*    right-hand side function in y'=f(t,y), the inital time T0, and */
+                         /*    the initial dependent variable vector y.  Note: since this */
+                         /*    problem is fully implicit, we set f_E to NULL and f_I to f. */
+                         flag = ARKodeInit(arkode_mem, NULL, f, T0, y);
+                         if (check_flag(&flag, "ARKodeInit", 1)) return 1;
+
                          return 0;
                        } |]
   putStrLn $ show res
diff --git a/packages/sundials/src/helpers.c b/packages/sundials/src/helpers.c
index 3498119..3e1a9a8 100644
--- a/packages/sundials/src/helpers.c
+++ b/packages/sundials/src/helpers.c
@@ -1,5 +1,12 @@
 #include <stdio.h>
-#include "helpers.h"
+#include <math.h>
+#include <arkode/arkode.h>                 /* prototypes for ARKODE fcts., consts. */
+#include <nvector/nvector_serial.h>        /* serial N_Vector types, fcts., macros */
+#include <sunmatrix/sunmatrix_dense.h>     /* access to dense SUNMatrix            */
+#include <sunlinsol/sunlinsol_dense.h>     /* access to dense SUNLinearSolver      */
+#include <arkode/arkode_direct.h>          /* access to ARKDls interface           */
+#include <sundials/sundials_types.h>       /* definition of type realtype          */
+#include <sundials/sundials_math.h>
 
 /* Check function return value...
     opt == 0 means SUNDIALS function allocates memory so check if
@@ -35,3 +42,30 @@ int check_flag(void *flagvalue, const char *funcname, int opt)
 
   return 0;
 }
+
+/* f routine to compute the ODE RHS function f(t,y). */
+int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
+{
+  realtype *rdata = (realtype *) user_data;   /* cast user_data to realtype */
+  realtype lamda = rdata[0];                  /* set shortcut for stiffness parameter */
+  realtype u = NV_Ith_S(y,0);                 /* access current solution value */
+
+  /* fill in the RHS function: "NV_Ith_S" accesses the 0th entry of ydot */
+  NV_Ith_S(ydot,0) = lamda*u + 1.0/(1.0+t*t) - lamda*atan(t);
+
+  return 0;                                   /* return with success */
+}
+
+/* Jacobian routine to compute J(t,y) = df/dy. */
+int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
+        void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
+{
+  realtype *rdata = (realtype *) user_data;   /* cast user_data to realtype */
+  realtype lamda = rdata[0];                  /* set shortcut for stiffness parameter */
+  realtype *Jdata = SUNDenseMatrix_Data(J);
+  
+  /* Fill in Jacobian of f: set the first entry of the data array to set the (0,0) entry */
+  Jdata[0] = lamda;
+
+  return 0;                                   /* return with success */
+}
diff --git a/packages/sundials/src/helpers.h b/packages/sundials/src/helpers.h
index d574ff5..ec11355 100644
--- a/packages/sundials/src/helpers.h
+++ b/packages/sundials/src/helpers.h
@@ -1 +1,27 @@
+#if defined(SUNDIALS_EXTENDED_PRECISION)
+#define GSYM "Lg"
+#define ESYM "Le"
+#define FSYM "Lf"
+#else
+#define GSYM "g"
+#define ESYM "e"
+#define FSYM "f"
+#endif
+
+/* Check function return value...
+    opt == 0 means SUNDIALS function allocates memory so check if
+             returned NULL pointer
+    opt == 1 means SUNDIALS function returns a flag so check if
+             flag >= 0
+    opt == 2 means function allocates memory so check if returned
+             NULL pointer  
+*/
 int check_flag(void *flagvalue, const char *funcname, int opt);
+
+/* f routine to compute the ODE RHS function f(t,y). */
+int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
+
+/* Jacobian routine to compute J(t,y) = df/dy. */
+int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
+        void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
+
diff --git a/packages/sundials/sundials.cabal b/packages/sundials/sundials.cabal
index 85a96fa..a68f33a 100644
--- a/packages/sundials/sundials.cabal
+++ b/packages/sundials/sundials.cabal
@@ -16,6 +16,8 @@ build-type:          Simple
 extra-source-files:  ChangeLog.md
 cabal-version:       >=1.10
 
+extra-source-files:  src/helpers.c, src/helpers.h
+
 executable sundials
   main-is:             Main.hs
   -- other-modules:       
@@ -28,3 +30,4 @@ executable sundials
   hs-source-dirs:      src
   default-language:    Haskell2010
   extra-libraries:     sundials_arkode
+  C-sources:           src/helpers.c src/helpers.h