Added cvode based solver. Still seems to have what I assume are numerical stability issues

src/Makefile

@@ -31,7 +31,8 @@ test1: app
 	#./grid2ode ../data/test3.txt
 	#./grid2ode ../data/test2.txt
 	./m2c circuit.mo
-	${CXX} ${CFLAGS} ${INCLUDE} -o sim func.cpp euler.cpp event.cpp 
+	#${CXX} ${CFLAGS} ${INCLUDE} -o sim func.cpp euler.cpp event.cpp 
+	${CXX} ${CFLAGS} ${INCLUDE} -o sim func.cpp cvode.cpp event.cpp -lsundials_cvode
 
 app: ${OBJS}
 	${CXX} -o grid2ode ${CFLAGS} ${OBJS} ${LIBS}

src/cvode.cpp

+#include "sim.h"
+#include <cfloat>
+#include <algorithm>
+#include <vector>
+#include <list>
+#include <cassert>
+#include <iostream>
+#include <cmath>
+#include <cstring>
+#include <cvode/cvode.h>
+#include <cvode/cvode_direct.h>
+#include <nvector/nvector_serial.h>
+#include <sunmatrix/sunmatrix_dense.h>
+#include <sunlinsol/sunlinsol_dense.h>
+#include <sundials/sundials_types.h>
+
+using namespace std;
+
+// Numerical tolerace. Picked these arbitrarily.
+static const double tolerance = 1E-10;
+static const double maxStep = 1E-8;
+// Model variable values
+static double *q, *dq, *alg;
+static double hmin = 0.0;
+
+static int f(realtype t, N_Vector y, N_Vector ydot, void* user_data)
+{
+	for (int i = 0; i < numStateVars; i++)
+		q[i] = NV_Ith_S(y,i);
+	event(dq,q,t);
+	func(dq,q,alg,t);
+	event(dq,q,t);
+	for (int i = 0; i < numStateVars; i++)
+		NV_Ith_S(ydot,i) = dq[i];
+	return 0;
+}
+
+void print(double t)
+{
+	cout << t;
+	for (int i = 0; i < numStateVars; i++)
+		cout << " " << q[i];
+	for (int i = 0; i < numAlgVars; i++)
+		cout << " " << alg[i];
+	cout << " " << hmin << endl;
+}
+
+int main(int argc, char** argv)
+{
+	double t = 0, cint = 0.0, tend = 0.0, tReport = 0.0;
+	N_Vector y, abstol;
+	SUNMatrix A;
+	SUNLinearSolver LS;
+	void *cvode_mem;
+
+	for (int i = 1; i < argc; i++)
+	{
+		if (strcmp("cint",argv[i]) == 0)
+			cint = atof(argv[++i]);
+		else if (strcmp("tend",argv[i]) == 0)
+			tend = atof(argv[++i]);
+	}
+	if (cint < 0.0 || tend <= 0.0)
+	{
+		cout << "Requires: cint <seconds> tend <seconds>" << endl;
+		return 0;
+	}
+	y = abstol = NULL;
+	A = NULL;
+	LS = NULL;
+	cvode_mem = NULL;
+
+	y = N_VNew_Serial(numStateVars);
+	abstol = N_VNew_Serial(numStateVars);
+
+	alg = new double[numAlgVars];
+	dq = new double[numStateVars];
+	q = new double[numStateVars];
+	// Calculate initial values
+	init_state(q);
+	for (int i = 0; i < numStateVars; i++)
+	{
+		NV_Ith_S(y,i) = q[i];
+		NV_Ith_S(abstol,i) = tolerance;
+	}
+	// Setup cvode	
+	cvode_mem = CVodeCreate(CV_BDF,CV_NEWTON);
+	CVodeSetMaxStep(cvode_mem,maxStep);
+	CVodeInit(cvode_mem,f,0.0,y);
+	CVodeSVtolerances(cvode_mem,tolerance,abstol);
+	A = SUNDenseMatrix(numStateVars,numStateVars);
+	LS = SUNDenseLinearSolver(y,A);
+	CVDlsSetLinearSolver(cvode_mem,LS,A);
+	// Report initial state
+	event(dq,q,t);
+	func(dq,q,alg,t);
+	event(dq,q,t);
+	print(t);
+	hmin = DBL_MAX;
+	// Run simulation
+	while (t < tend)
+	{
+		double tPrev = t;
+		CVode(cvode_mem,t+maxStep,y,&t,CV_ONE_STEP);
+		hmin = ::min(hmin,t-tPrev);
+		if (t >= tReport+cint)
+		{
+			tReport += cint;
+			for (int i = 0; i < numStateVars; i++)
+				q[i] = NV_Ith_S(y,i);
+			event(dq,q,t);
+			func(dq,q,alg,t);
+			event(dq,q,t);
+			print(t);
+			hmin = DBL_MAX;
+		}
+	}
+	// Cleanup and exit
+	delete [] dq;
+	delete [] q;
+	delete [] alg;
+
+	N_VDestroy(y);
+	N_VDestroy(abstol);
+	CVodeFree(&cvode_mem);
+	SUNLinSolFree(LS);
+	SUNMatDestroy(A);
+	return 0;
+}