diff --git a/examples/oregonator/SConscript b/examples/oregonator/SConscript
new file mode 100644
index 0000000..f977215
--- /dev/null
+++ b/examples/oregonator/SConscript
@@ -0,0 +1,23 @@
+Import('env')
+cObj = []
+c_src = Glob('*.c')
+import os
+
+FD = env['FINITE_DIFFERENCE']
+
+for src in c_src:
+    if FD and 'jacob' in str(src):
+        continue
+    cObj.append(env.Object(src, variant_dir=env['variant']))
+
+
+cudaObj = []
+if env['build_cuda']:
+    cuda_src = Glob('*.cu')
+
+    for src in cuda_src:
+        if FD and 'jacob' in str(src):
+            continue
+        cudaObj.append(env.CUDAObject(src, variant_dir=env['variant']))
+
+Return ('cObj', 'cudaObj')
\ No newline at end of file
diff --git a/examples/oregonator/dydt.c b/examples/oregonator/dydt.c
new file mode 100644
index 0000000..93c2fa4
--- /dev/null
+++ b/examples/oregonator/dydt.c
@@ -0,0 +1,43 @@
+/**
+ * \file
+ * \brief An implementation of the Oregonator right hand side (y' = f(y)) function.
+ *
+ * Implements the evaluation of the derivative of the state vector with respect to time, \f$\dot{\vec{y}}\f$
+ *
+ */
+
+#include "header.h"
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator {
+#endif
+
+/**
+ * \brief An implementation of the RHS of the Oregonator equation
+ * \param[in]        t         The current system time
+ * \param[in]        mu        Dummy parameter, needed for compatibility only
+ * \param[in]        y         The state vector
+ * \param[out]       dy        The output RHS (dydt) vector
+ *
+ * The `y` and `dy` vectors supplied here are local versions of the global state vectors.
+ * They have been transformed from the global Column major (Fortan) ordering to a local 1-D vector
+ * Hence the vector accesses can be done in a simple manner below, i.e. y[0] -> \f$y_1\f$, y[1] -> \f$y_2\f$, etc.
+ * @see solver_generic.c
+ */
+void dydt (const double t, const double mu, const double * __restrict__ y, double * __restrict__ dy) {
+
+  double s = 77.27;
+  double q = 8.375E-6;
+  double w = 0.161;
+
+  dy[0] = s * (y[0] - y[0] * y[1] + y[1] - (q * y[0]* y[0]));
+  dy[1] = (y[2] - y[1] - y[0] * y[1]) / s;
+  dy[2] = w * (y[0] - y[2]);
+
+} // end dydt
+
+
+#ifdef GENERATE_DOCS
+}
+#endif
diff --git a/examples/oregonator/dydt.cu b/examples/oregonator/dydt.cu
new file mode 100644
index 0000000..41838fc
--- /dev/null
+++ b/examples/oregonator/dydt.cu
@@ -0,0 +1,47 @@
+/**
+ * \file
+ * \brief A CUDA implementation of the Oregonator right hand side (y' = f(y)) function.
+ *
+ * Implements the CUDA evaluation of the derivative of the state vector with respect to time, \f$\dot{\vec{y}}\f$
+ *
+ */
+
+#include "header.cuh"
+#include "gpu_macros.cuh"
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator_cu {
+#endif
+
+/**
+ * \brief An implementation of the RHS of the Oregonator equation
+ * \param[in]        t         The current system time
+ * \param[in]        mu        Dummy parameter, needed for compatibility only
+ * \param[in]        y         The state vector
+ * \param[out]       dy        The output RHS (dydt) vector
+ * \param[in]        d_mem     The mechanism_memory struct.  In future versions, this will be used to access the \f$\mu\f$ parameter to have a consistent interface.
+ *
+ * The `y` and `dy` vectors supplied here are the global state vectors.
+ * Hence the vector accesses must be done with the global thread ID.
+ * The gpu_macros.cuh file defines some useful macros to simplify indexing
+ * @see gpu_macros.cuh
+ */
+ __device__
+void dydt (const double t, const double mu, const double * __restrict__ y, double * __restrict__ dy,
+           const mechanism_memory * __restrict__ d_mem) {
+
+   double s = 77.27;
+   double q = 8.375E-6;
+   double w = 0.161;
+
+   dy[INDEX(0)] = s * (y[INDEX(0)] - y[INDEX(0)] * y[INDEX(1)] + y[INDEX(1)] - (q * y[INDEX(0)]* y[INDEX(0)]));
+   dy[INDEX(1)] = (y[INDEX(2)] - y[INDEX(1)] - y[INDEX(0)] * y[INDEX(1)]) / s;
+   dy[INDEX(2)] = w * (y[INDEX(0)] - y[INDEX(2)]);
+
+} // end dydt
+
+
+#ifdef GENERATE_DOCS
+}
+#endif
diff --git a/examples/oregonator/dydt.cuh b/examples/oregonator/dydt.cuh
new file mode 100644
index 0000000..8785d1c
--- /dev/null
+++ b/examples/oregonator/dydt.cuh
@@ -0,0 +1,34 @@
+/**
+ * \file
+ * \brief Contains header definitions for the CUDA RHS function for the Oregonator example
+ *
+ */
+
+#ifndef DYDT_CUH
+#define DYDT_CUH
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator_cu {
+#endif
+
+
+
+/**
+ * \brief An implementation of the RHS of the Oregonator equation
+ * \param[in]        t         The current system time
+ * \param[in]        mu        Dummy parameter, needed for compatibility only
+ * \param[in]        y         The state vector
+ * \param[out]       dy        The output RHS (dydt) vector
+ * \param[in]        d_mem     The mechanism_memory struct.  In future versions, this will be used to access the \f$\mu\f$ parameter to have a consistent interface.
+ *
+ */
+__device__
+void dydt (const double t, const double mu, const double * __restrict__ y, double * __restrict__ dy,
+           const mechanism_memory * __restrict__ d_mem);
+
+#ifdef GENERATE_DOCS
+}
+#endif
+
+#endif
diff --git a/examples/oregonator/dydt.h b/examples/oregonator/dydt.h
new file mode 100644
index 0000000..b9eae10
--- /dev/null
+++ b/examples/oregonator/dydt.h
@@ -0,0 +1,30 @@
+/**
+ * \file
+ * \brief Contains header definitions for the RHS function for the Oregonator example
+ *
+ */
+
+#ifndef DYDT_H
+#define DYDT_H
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator {
+#endif
+
+
+
+/**
+ * \brief An implementation of the RHS of the Oregonator equation
+ * \param[in]        t         The current system time
+ * \param[in]        mu        Dummy parameter, needed for compatibility only
+ * \param[in]        y         The state vector
+ * \param[out]       dy        The output RHS (dydt) vector
+ */
+void dydt (const double t, const double mu, const double * __restrict__ y, double * __restrict__ dy);
+
+#ifdef GENERATE_DOCS
+}
+#endif
+
+#endif
diff --git a/examples/oregonator/gpu_macros.cuh b/examples/oregonator/gpu_macros.cuh
new file mode 100644
index 0000000..0bd4fac
--- /dev/null
+++ b/examples/oregonator/gpu_macros.cuh
@@ -0,0 +1,41 @@
+/**
+ * \file
+ * \brief Defines some simple macros to simplify GPU indexing
+ *
+ */
+
+#ifndef GPU_MACROS_CUH
+#define GPU_MACROS_CUH
+#include <stdio.h>
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include <helper_cuda.h>
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator_cu {
+#endif
+
+//! The total number of threads in the Grid, provides an offset between vector entries
+#define GRID_DIM (blockDim.x * gridDim.x)
+//! The global CUDA thread index
+#define T_ID (threadIdx.x + blockIdx.x * blockDim.x)
+//! Convenience macro to get the value of a vector at index i, calculated as i * #GRID_DIM + #T_ID
+#define INDEX(i) (T_ID + (i) * GRID_DIM)
+
+#define cudaErrorCheck(ans) { gpuAssert((ans), __FILE__, __LINE__); }
+inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=true)
+{
+    if (code != cudaSuccess)
+    {
+        fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
+        if (abort) exit(code);
+    }
+}
+
+
+#ifdef GENERATE_DOCS
+}
+#endif
+
+#endif
diff --git a/examples/oregonator/gpu_memory.cu b/examples/oregonator/gpu_memory.cu
new file mode 100644
index 0000000..7cab7d6
--- /dev/null
+++ b/examples/oregonator/gpu_memory.cu
@@ -0,0 +1,68 @@
+/*!
+ * \file
+ * \brief Initializes and calculates required GPU memory
+ */
+
+#include "gpu_memory.cuh"
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator_cu {
+#endif
+
+/**
+ * \brief Calculates and returns the total memory size (in bytes) required by an individual thread for the
+ *        mechanism_memory struct.
+ */
+size_t required_mechanism_size() {
+  //returns the total required size for the mechanism per thread
+  size_t mech_size = 0;
+  //state vector y
+  mech_size += NSP;
+  //dydt vector
+  mech_size += NSP;
+  //Jacobian
+  mech_size += NSP * NSP;
+  //and mu parameter
+  mech_size += 1;
+  return mech_size * sizeof(double);
+}
+/**
+ * \brief Initializes the host and device mechanism_memory structs.  This is required in order to enable
+ *        passing the struct to CUDA
+ * \param[in]               padded              The padded number of threads to be used by the CUDA solver
+ * \param[in,out]           h_mem               The host version of the mechanism_memory struct to initialize.
+ * \param[in,out]           d_mem               The device version of the mechanism_memory struct to copy the resulting host mechanism_memory struct to.
+ */
+void initialize_gpu_memory(int padded, mechanism_memory** h_mem, mechanism_memory** d_mem)
+{
+  // Allocate storage for the device struct
+  cudaErrorCheck( cudaMalloc(d_mem, sizeof(mechanism_memory)) );
+  //allocate the device arrays on the host pointer
+  cudaErrorCheck( cudaMalloc(&((*h_mem)->y), NSP * padded * sizeof(double)) );
+  cudaErrorCheck( cudaMalloc(&((*h_mem)->dy), NSP * padded * sizeof(double)) );
+  cudaErrorCheck( cudaMalloc(&((*h_mem)->var), 1 * padded * sizeof(double)) );
+  cudaErrorCheck( cudaMalloc(&((*h_mem)->jac), NSP * NSP * padded * sizeof(double)) );
+  // set non-initialized values to zero
+  cudaErrorCheck( cudaMemset((*h_mem)->dy, 0, NSP * padded * sizeof(double)) );
+  cudaErrorCheck( cudaMemset((*h_mem)->jac, 0, NSP * NSP * padded * sizeof(double)) );
+  // and copy to device pointer
+  cudaErrorCheck( cudaMemcpy(*d_mem, *h_mem, sizeof(mechanism_memory), cudaMemcpyHostToDevice) );
+}
+/**
+ * \brief Frees the host and device mechanism_memory structs
+ * \param[in,out]           h_mem               The host version of the mechanism_memory struct.
+ * \param[in,out]           d_mem               The device version of the mechanism_memory struct.
+ */
+void free_gpu_memory(mechanism_memory** h_mem, mechanism_memory** d_mem)
+{
+  cudaErrorCheck(cudaFree((*h_mem)->y));
+  cudaErrorCheck(cudaFree((*h_mem)->dy));
+  cudaErrorCheck(cudaFree((*h_mem)->var));
+  cudaErrorCheck(cudaFree((*h_mem)->jac));
+  cudaErrorCheck(cudaFree(*d_mem));
+}
+
+#ifdef GENERATE_DOCS
+}
+#endif
diff --git a/examples/oregonator/gpu_memory.cuh b/examples/oregonator/gpu_memory.cuh
new file mode 100644
index 0000000..edac909
--- /dev/null
+++ b/examples/oregonator/gpu_memory.cuh
@@ -0,0 +1,66 @@
+/*!
+ * \file
+ * \brief Headers for GPU memory initialization
+ */
+
+<<<<<<< HEAD
+=======
+#ifndef GPU_MEM_CUH
+#define GPU_MEM_CUH
+
+>>>>>>> master
+#include "header.cuh"
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator_cu {
+#endif
+
+/**
+ * \brief This struct is used to store memory for the CUDA RHS and Jacobian evaluation.
+ *        Along with the solver_memory struct, this must be initialized and passed to the solvers
+ *        the run.  @see solver_memory, initialize_gpu_memory, required_mechanism_size
+ *
+ * \param           y               The global state vector arrays
+ * \param           dy              The global dydt vector arrays
+ * \param           var             The global param array [Used for \f$\mu\f$]
+ * \param           jac             The global Jacobian arrays
+ *
+ * This is also heavily used when using with @pyJac to hold additional arrays
+ * for evaluating chemical kinetic reactions.
+ */
+struct mechanism_memory {
+  double * y;
+  double * dy;
+  double * var;
+  double * jac;
+};
+
+/**
+ * \brief Calculates and returns the total memory size (in bytes) required by an individual thread for the
+ *        mechanism_memory struct.
+ */
+size_t required_mechanism_size();
+/**
+ * \brief Initializes the host and device mechanism_memory structs.  This is required in order to enable
+ *        passing the struct to CUDA
+ * \param[in]               padded              The padded number of threads to be used by the CUDA solver
+ * \param[in,out]           h_mem               The host version of the mechanism_memory struct to initialize.
+ * \param[in,out]           d_mem               The device version of the mechanism_memory struct to copy the resulting host mechanism_memory struct to.
+ */
+void initialize_gpu_memory(int padded, mechanism_memory** h_mem, mechanism_memory** d_mem);
+/**
+ * \brief Frees the host and device mechanism_memory structs
+ * \param[in,out]           h_mem               The host version of the mechanism_memory struct.
+ * \param[in,out]           d_mem               The device version of the mechanism_memory struct.
+ */
+void free_gpu_memory(mechanism_memory** h_mem, mechanism_memory** d_mem);
+
+#ifdef GENERATE_DOCS
+}
+#endif
+<<<<<<< HEAD
+=======
+
+#endif
+>>>>>>> master
diff --git a/examples/oregonator/header.cuh b/examples/oregonator/header.cuh
new file mode 100644
index 0000000..6189255
--- /dev/null
+++ b/examples/oregonator/header.cuh
@@ -0,0 +1,56 @@
+/*!
+ * \file
+ * \brief An example header file that defines system size, memory functions and other required methods
+ *        for integration of the Oregonator's equation with the CUDA solvers
+ */
+
+#ifndef HEADER_GUARD_CUH
+#define HEADER_GUARD_CUH
+
+#include <stdlib.h>
+#include "gpu_macros.cuh"
+#include "gou_memory.cuh"
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator_cu {
+#endif
+
+//! The IVP system size
+#define NSP (3)
+//! Input vector size (in read_initial_conditions)
+#define NN (NSP)
+
+/*!
+ *
+ * \brief Set same ICs for all problems
+ * \param NUM       The number of initial value problems
+ * \param y_host    The state vectors to initialize
+ * \param var_host  The vector of \f$mu\f$ parameters for the Oregonator equation, not working now
+ *
+ */
+void set_same_initial_conditions(int NUM, double** y_host, double** var_host);
+
+//dummy definitions that are used for pyJac
+
+/**
+ * \brief Not needed for Oregonator
+ *
+ *  In pyJac, these are used to transform the input/output vectors to deal with moving the
+ *         last species mass fraction
+ */
+void apply_mask(double* y_host);
+/**
+ * \brief Not needed for Oregonator
+ *
+ *  In pyJac, these are used to transform the input/output vectors to deal with moving the
+ *         last species mass fraction
+ */
+void apply_reverse_mask(double* y_host);
+
+
+#ifdef GENERATE_DOCS
+}
+#endif
+
+#endif
diff --git a/examples/oregonator/header.h b/examples/oregonator/header.h
new file mode 100644
index 0000000..737388a
--- /dev/null
+++ b/examples/oregonator/header.h
@@ -0,0 +1,62 @@
+/*!
+ * \file
+ * \brief An example header file that defines system size and other required methods
+ *        for integration of the Oregonator's equation.
+ */
+
+#ifndef HEADER_GUARD_H
+#define HEADER_GUARD_H
+
+#include <stdlib.h>
+#include <math.h>
+
+//include OpenMP if available
+#ifdef _OPENMP
+ #include <omp.h>
+#else
+ #define omp_get_max_threads() 1
+ #define omp_get_num_threads() 1
+#endif
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator {
+#endif
+
+//! The IVP system size
+#define NSP (3)
+//! Input vector size (in read_initial_conditions)
+#define NN (NSP)
+
+/*!
+ *
+ * \brief Set same ICs for all problems
+ * \param NUM       The number of initial value problems
+ * \param y_host    The state vectors to initialize
+ * \param var_host  The vector of \f$mu\f$ parameters for the Oregonator equation, not currently working
+ *
+ */
+void set_same_initial_conditions(int NUM, double** y_host, double** var_host);
+
+//dummy definitions that are used for pyJac
+
+/**
+ * \brief Not needed for Oregonator
+ *
+ *  In pyJac, these are used to transform the input/output vectors to deal with moving the
+ *         last species mass fraction
+ */
+void apply_mask(double* y_host);
+/**
+ * \brief Not needed for Oregonator
+ *
+ *  In pyJac, these are used to transform the input/output vectors to deal with moving the
+ *         last species mass fraction
+ */
+void apply_reverse_mask(double* y_host);
+
+#ifdef GENERATE_DOCS
+}
+#endif
+
+#endif
diff --git a/examples/oregonator/ics.c b/examples/oregonator/ics.c
new file mode 100644
index 0000000..26a86dc
--- /dev/null
+++ b/examples/oregonator/ics.c
@@ -0,0 +1,58 @@
+/**
+ * \file
+ * \brief Sets same Initial Conditions (ICs) for all problems
+ *
+ * This provides simple definitions of the set_same_initial_conditions function
+ * used when SAME_ICS is defined
+ */
+
+#include "header.h"
+
+#ifdef GENERATE_DOCS
+//put this in the van der Pol namespace for documentation
+namespace oregonator {
+#endif
+
+
+/*!
+ *
+ * \brief Set same ICs for all problems
+ * \param NUM       The number of initial value problems
+ * \param y_host    The state vectors to initialize
+ * \param var_host  The vector of \f$\mu\f$ parameters for the Oregonator equation, not currently working
+ *
+ */
+void set_same_initial_conditions(int NUM, double** y_host, double** var_host)
+{
+    //init vectors
+    (*y_host) = (double*)malloc(NUM * NSP * sizeof(double));
+    (*var_host) = (double*)malloc(NUM * sizeof(double));
+    //now set the values
+    for (int i = 0; i < NUM; ++i){
+      //Dummy value, this won't do anything important for the Oregonator
+      (*var_host)[i] = 1000;
+      //set y values
+      (*y_host)[i] = 1;
+      (*y_host)[i + NUM] = 1;
+      (*y_host)[i + 2 * NUM] = 2;
+    }
+}
+
+/**
+ * \brief Not needed for Oregonator
+ *
+ *  In pyJac, these are used to transform the input/output vectors to deal with moving the
+ *         last species mass fraction
+ */
+void apply_mask(double* y_host) {}
+/**
+ * \brief Not needed for Oregonator
+ *
+ *  In pyJac, these are used to transform the input/output vectors to deal with moving the
+ *         last species mass fraction
+ */
+void apply_reverse_mask(double* y_host) {}
+
+#ifdef GENERATE_DOCS
+}
+#endif
diff --git a/examples/oregonator/ics.cu b/examples/oregonator/ics.cu
new file mode 100644
index 0000000..10e0310
--- /dev/null
+++ b/examples/oregonator/ics.cu
@@ -0,0 +1,58 @@
+/**
+ * \file
+ * \brief Sets same Initial Conditions (ICs) for all problems
+ *
+ * This provides simple definitions of the set_same_initial_conditions function
+ * used when SAME_ICS is defined
+ */
+
+#include "header.h"
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator_cu {
+#endif
+
+
+/*!
+ *
+ * \brief Set same ICs for all problems
+ * \param NUM       The number of initial value problems
+ * \param y_host    The state vectors to initialize
+ * \param var_host  The vector of \f$\mu\f$ parameters for the Oregonator equation
+ *
+ */
+void set_same_initial_conditions(int NUM, double** y_host, double** var_host)
+{
+    //init vectors
+    (*y_host) = (double*)malloc(NUM * NSP * sizeof(double));
+    (*var_host) = (double*)malloc(NUM * sizeof(double));
+    //now set the values
+    for (int i = 0; i < NUM; ++i){
+      //Dummy value, this won't do anything important for the Oregonator
+      (*var_host)[i] = 1000;
+      //set y values
+      (*y_host)[i] = 1;
+      (*y_host)[i + NUM] = 1;
+      (*y_host)[i + 2 * NUM] = 2;
+    }
+}
+
+/**
+ * \brief Not needed for Oregonator
+ *
+ *  In pyJac, these are used to transform the input/output vectors to deal with moving the
+ *         last species mass fraction
+ */
+void apply_mask(double* y_host) {}
+/**
+ * \brief Not needed for Oregonator
+ *
+ *  In pyJac, these are used to transform the input/output vectors to deal with moving the
+ *         last species mass fraction
+ */
+void apply_reverse_mask(double* y_host) {}
+
+#ifdef GENERATE_DOCS
+}
+#endif
diff --git a/examples/oregonator/jacob.c b/examples/oregonator/jacob.c
new file mode 100644
index 0000000..96a66ad
--- /dev/null
+++ b/examples/oregonator/jacob.c
@@ -0,0 +1,82 @@
+/**
+ * \file
+ * \brief An implementation of the Oregonator jacobian \f$\frac{\partial \dot{\vec{y}}}{\partial \vec{y}}\f$
+ *
+ * Implementes the evaluation of Oregonator Jacobian
+ *
+ */
+
+#include "header.h"
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator {
+#endif
+
+#define TRANSPOSE
+
+/**
+ * \brief An implementation of the Oregonator jacobian
+ *
+ * \param[in]           t               The current system time
+ * \param[in]           mu              Dummy parameter, needed for compatibility only
+ * \param[in]           y               The state vector at time t
+ * \param[out]          jac             The jacobian to populate
+ *
+ *  The Jacobian is in a local Column-major (Fortran) order.  As with dydt(), this function operates on local
+ *  copies of the global state vector and jacobian.  Hence simple linear indexing can be used here.
+ *  @see solver_generic.c
+ *
+ */
+void eval_jacob (const double t, const double mu, const double * __restrict__ y, double * __restrict__ jac)
+{
+
+  double s = 77.27;
+  double q = 8.375E-6;
+  double w = 0.161;
+
+#ifndef TRANSPOSE
+  // dydot_0 / dy0
+  jac[0] = s * (-1 * y[1] + 1 - q * 2 * y[0]);
+  // dydot_0 / dy1
+  jac[1] = s * (1 - y[0]);
+  // dydot_0 / dy2
+  jac[2] = 0;
+  // dydot_1 / dy0
+  jac[3] = -1 * y[1] / s;
+  // dydot_1 / dy1
+  jac[4] = (-1 - y[0]) / s;
+  // dydot_1 / dy2
+  jac[5] = 1 / s;
+  // dydot_2 / dy0
+  jac[6] = w;
+  // dydot_2 / dy1
+  jac[7] = 0;
+  // dydot_2 / dy2
+  jac[8] = -1 * w;
+#else
+  // dydot_0 / dy0
+  jac[0] = s * (-1 * y[1] + 1 - q * 2 * y[0]);
+  // dydot_1 / dy0
+  jac[1] = -1 * y[1] / s;
+  // dydot_2 / dy0
+  jac[2] = w;
+  // dydot_0 / dy1
+  jac[3] = s * (1 - y[0]);
+  // dydot_1 / dy1
+  jac[4] = (-1 - y[0]) / s;
+  // dydot_2 / dy1
+  jac[5] = 0;
+  // dydot_0 / dy2
+  jac[6] = 0;
+  // dydot_1 / dy2
+  jac[7] = 1 / s;
+  // dydot_2 / dy2
+  jac[8] = -1 * w;
+#endif
+
+}
+
+#ifdef GENERATE_DOCS
+}
+#endif
diff --git a/examples/oregonator/jacob.cu b/examples/oregonator/jacob.cu
new file mode 100644
index 0000000..7c55b13
--- /dev/null
+++ b/examples/oregonator/jacob.cu
@@ -0,0 +1,47 @@
+/**
+ * \file
+ * \brief A CUDA implementation of the Oregonator jacobian \f$\frac{\partial \dot{\vec{y}}}{\partial \vec{y}}\f$
+ *
+ * Implementes the CUDA evaluation of Oregonator Jacobian
+ *
+ */
+
+#include "header.cuh"
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator_cu {
+#endif
+
+/**
+ * \brief An implementation of the Oregonator jacobian
+ *
+ * \param[in]           t               The current system time
+ * \param[in]           mu              Dummy parameter, needed for compatibility only
+ * \param[in]           y               The state vector at time t
+ * \param[out]          jac             The jacobian to populate
+ * \param[in]           d_mem           The mechanism_memory struct.  In future versions, this will be used to access the \f$\mu\f$ parameter to have a consistent interface.
+ *
+ *  The Jacobian is in Column-major (Fortran) order.  As with dydt(), this function operates directly on
+ *  global state vector and jacobian.  Hence we use the #INDEX macro defined in gpu_macros.cuh here.
+ *  @see dydt()
+ *
+ */
+__device__
+void eval_jacob (const double t, const double mu, const double * __restrict__ y, double * __restrict__ jac, const mechanism_memory * __restrict__ d_mem)
+{
+  jac[INDEX(0)] = s * (-1 * y[INDEX(1)] + 1 - q * 2 * y[INDEX(0)]);
+  jac[INDEX(1)] = s * (1 - y[INDEX(0)]);
+  jac[INDEX(2)] = 0;
+  jac[INDEX(3)] = -1 * y[INDEX(1)] / s;
+  jac[INDEX(4)] = (-1 - y[INDEX(0)]) / s;
+  jac[INDEX(5)] = 1 / s;
+  jac[INDEX(6)] = w;
+  jac[INDEX(7)] = 0;
+  jac[INDEX(8)] = -1 * w;
+
+}
+
+#ifdef GENERATE_DOCS
+}
+#endif
diff --git a/examples/oregonator/jacob.cuh b/examples/oregonator/jacob.cuh
new file mode 100644
index 0000000..0b642aa
--- /dev/null
+++ b/examples/oregonator/jacob.cuh
@@ -0,0 +1,28 @@
+/**
+ * \file
+ * \brief Contains a header definition for the CUDA van der Pol Jacobian evaluation
+ *
+ */
+
+#ifdef GENERATE_DOCS
+//put this in the van der Pol namespace for documentation
+namespace van_der_pol_cu {
+#endif
+
+
+/**
+ * \brief An implementation of the van der Pol jacobian
+ *
+ * \param[in]           t               The current system time
+ * \param[in]           mu              The van der Pol parameter
+ * \param[in]           y               The state vector at time t
+ * \param[out]          jac             The jacobian to populate
+ * \param[in]           d_mem           The mechanism_memory struct.  In future versions, this will be used to access the \f$\mu\f$ parameter to have a consistent interface.
+ *
+ */
+__device__
+void eval_jacob (const double t, const double mu, const double * __restrict__ y, double * __restrict__ jac, const mechanism_memory * __restrict__ d_mem);
+
+#ifdef GENERATE_DOCS
+}
+#endif
\ No newline at end of file
diff --git a/examples/oregonator/jacob.h b/examples/oregonator/jacob.h
new file mode 100644
index 0000000..9fcd8d1
--- /dev/null
+++ b/examples/oregonator/jacob.h
@@ -0,0 +1,32 @@
+/**
+ * \file
+ * \brief Contains a header definition for the Oregonator Jacobian evaluation
+ *
+ */
+
+#ifndef JACOB_H
+#define JACOB_H
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator {
+#endif
+
+
+/**
+ * \brief An implementation of the Oregonator jacobian
+ *
+ * \param[in]           t               The current system time
+ * \param[in]           mu              Dummy parameter, needed for compatibility only
+ * \param[in]           y               The state vector at time t
+ * \param[out]          jac             The jacobian to populate
+ *
+ */
+void eval_jacob (const double t, const double mu, const double * __restrict__ y, double * __restrict__ jac);
+
+#ifdef GENERATE_DOCS
+}
+
+#endif
+
+#endif
diff --git a/examples/oregonator/launch_bounds.cuh b/examples/oregonator/launch_bounds.cuh
new file mode 100644
index 0000000..cce4de4
--- /dev/null
+++ b/examples/oregonator/launch_bounds.cuh
@@ -0,0 +1,16 @@
+/**
+ *  \file
+ *  \brief A number of definitions that control CUDA kernel launches
+ */
+
+#ifndef LAUNCH_BOUNDS_CUH
+#define LAUNCH_BOUNDS_CUH
+
+//! The target number of threads per block
+#define TARGET_BLOCK_SIZE (64)
+//! Shared memory per thread
+#define SHARED_SIZE (0)
+//! Large L1 cache active
+#define PREFERL1
+
+#endif
\ No newline at end of file
diff --git a/examples/oregonator/oregonator.dox b/examples/oregonator/oregonator.dox
new file mode 100644
index 0000000..0159883
--- /dev/null
+++ b/examples/oregonator/oregonator.dox
@@ -0,0 +1,5 @@
+/**
+\page oregonator The Oregonator problem
+
+More to come!
+*/
diff --git a/examples/oregonator/plotter.py b/examples/oregonator/plotter.py
new file mode 100644
index 0000000..71ccdda
--- /dev/null
+++ b/examples/oregonator/plotter.py
@@ -0,0 +1,39 @@
+# imports
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+
+# load log files
+files = [os.path.join('log', file) for file in os.listdir('./log/')
+         if file.endswith('.bin')]
+files = [file for file in files if os.path.isfile(file)]
+
+linestyles = ['-', '--', '-.']
+colorwheel = ['r', 'g', 'b', 'k']
+
+# load data
+for i, file in enumerate(files):
+    arr = np.fromfile(file)
+    # reshape to a 2001 x 3 matrix (time, y1, y2)
+    arr = np.reshape(arr, (2001, 3))
+
+    # get the solver name
+    label = file.split('-')[0]
+    label = label[label.index('/') + 1:]
+    label += '-GPU' if 'gpu' in file else ''
+
+    # plot y1
+    plt.plot(arr[:, 0][::10], arr[:, 1][::10],
+             linestyle=linestyles[i % len(linestyles)],
+             label=label, color=colorwheel[i % len(colorwheel)])
+
+# make legend
+plt.legend(loc=0)
+
+# title and labels
+plt.title('Oregonator equation')
+plt.xlabel('Time(s)')
+plt.ylabel('$y_1$')
+
+# and save fig
+plt.savefig('oregonator.png', dpi=300, size=(5, 3))
diff --git a/examples/oregonator/sparse_multiplier.c b/examples/oregonator/sparse_multiplier.c
new file mode 100644
index 0000000..864956c
--- /dev/null
+++ b/examples/oregonator/sparse_multiplier.c
@@ -0,0 +1,31 @@
+/**
+ * \file
+ * \brief Implementation for Jacobian vector multiplication, used in exponential integrators
+ *
+ */
+
+#include "sparse_multiplier.h"
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator {
+#endif
+
+
+/**
+ * \brief Implements Jacobian \ vector multiplication in sparse (or unrolled) form
+ * \param[in]           A           The (NSP x NSP) Jacobian matrix, see eval_jacob() for details on layout
+ * \param[in]           Vm          The (NSP x 1) vector to multiply by
+ * \param[out]          w           The (NSP x 1) vector to store the result in, \f$w := A * Vm\f$
+ */
+ //Needs to be modified for Oregonator
+void sparse_multiplier(const double * A, const double * Vm, double* w) {
+  w[0] =  A[0] * Vm[0] +  A[3] * Vm[1] + A[6] * Vm[2];
+  w[1] =  A[1] * Vm[0] +  A[4] * Vm[1] + A[7] * Vm[2];
+  w[2] =  A[2] * Vm[0] +  A[5] * Vm[1] + A[8] * Vm[2];
+}
+
+
+#ifdef GENERATE_DOCS
+}
+#endif
diff --git a/examples/oregonator/sparse_multiplier.cu b/examples/oregonator/sparse_multiplier.cu
new file mode 100644
index 0000000..37cad4d
--- /dev/null
+++ b/examples/oregonator/sparse_multiplier.cu
@@ -0,0 +1,31 @@
+/**
+ * \file
+ * \brief Implementation for CUDA Jacobian vector multiplication, used in exponential integrators
+ *
+ */
+
+#include "sparse_multiplier.cuh"
+
+#ifdef GENERATE_DOCS
+//put this in the Oregonator namespace for documentation
+namespace oregonator_cu {
+#endif
+
+
+/**
+ * \brief Implements Jacobian \ vector multiplication in sparse (or unrolled) form
+ * \param[in]           A           The (NSP x NSP) Jacobian matrix, see eval_jacob() for details on layout
+ * \param[in]           Vm          The (NSP x 1) vector to multiply by
+ * \param[out]          w           The (NSP x 1) vector to store the result in, \f$w := A * Vm\f$
+ */
+__device__
+void sparse_multiplier(const double * A, const double * Vm, double* w) {
+  w[INDEX(0)] =  A[INDEX(0)] * Vm[INDEX(0)] +  A[INDEX(3)] * Vm[INDEX(1)] + A[INDEX(6)] * Vm[INDEX(2)];
+  w[INDEX(1)] =  A[INDEX(1)] * Vm[INDEX(0)] +  A[INDEX(4)] * Vm[INDEX(1)] + A[INDEX(7)] * Vm[INDEX(2)];
+  w[INDEX(2)] =  A[INDEX(2)] * Vm[INDEX(0)] +  A[INDEX(5)] * Vm[INDEX(1)] + A[INDEX(8)] * Vm[INDEX(2)];
+}
+
+
+#ifdef GENERATE_DOCS
+}
+#endif
diff --git a/examples/oregonator/sparse_multiplier.cuh b/examples/oregonator/sparse_multiplier.cuh
new file mode 100644
index 0000000..71181b3
--- /dev/null
+++ b/examples/oregonator/sparse_multiplier.cuh
@@ -0,0 +1,32 @@
+/**
+ * \file
+ * \brief Header definition for CUDA Jacobian vector multiplier, used in exponential integrators
+ *
+ */
+
+#ifndef SPARSE_HEAD_CUH
+#define SPARSE_HEAD_CUH
+
+#include "header.cuh"
+
+
+#ifdef GENERATE_DOCS
+//put this in the van der Pol namespace for documentation
+namespace van_der_pol_cu {
+#endif
+
+/**
+ * \brief Implements Jacobian \ vector multiplication in sparse (or unrolled) form
+ * \param[in]           A           The (NSP x NSP) Jacobian matrix, see eval_jacob() for details on layout
+ * \param[in]           Vm          The (NSP x 1) vector to multiply by
+ * \param[out]          w           The (NSP x 1) vector to store the result in, \f$w := A * Vm\f$
+ */
+__device__
+void sparse_multiplier (const double * A, const double * Vm, double* w);
+
+
+#ifdef GENERATE_DOCS
+}
+#endif
+
+#endif
diff --git a/examples/oregonator/sparse_multiplier.h b/examples/oregonator/sparse_multiplier.h
new file mode 100644
index 0000000..1db49a7
--- /dev/null
+++ b/examples/oregonator/sparse_multiplier.h
@@ -0,0 +1,25 @@
+/**
+ * \file
+ * \brief Header definition for Jacobian vector multiplier, used in exponential integrators
+ *
+ */
+
+#ifndef SPARSE_HEAD
+#define SPARSE_HEAD
+
+#include "header.h"
+
+
+#ifdef GENERATE_DOCS
+//put this in the van der Pol namespace for documentation
+namespace van_der_pol {
+#endif
+
+void sparse_multiplier (const double *, const double *, double*);
+
+
+#ifdef GENERATE_DOCS
+}
+#endif
+
+#endif