diff --git a/data/maref.lisp b/data/maref.lisp
index a7f1d6b3b1ea37c12ca04023fd99419648c4c7d4..deef9f636e6417bfed7a818e3159fbeb32e9c1b9 100644
--- a/data/maref.lisp
+++ b/data/maref.lisp
@@ -1,6 +1,6 @@
 ;; Get/set array or elements: cl-array, maref
 ;; Liam Healy 2008-08-27 22:43:10EDT maref.lisp
-;; Time-stamp: <2009-01-12 09:43:16EST maref.lisp>
+;; Time-stamp: <2009-01-17 17:15:19EST maref.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -89,8 +89,8 @@
 ;;; callback functions in solve-minimize-fit and in #'cl-array method
 ;;; for pointers.
 
-(eval-when (:compile-toplevel)
-(defun maref-function-picker (type-symbol category ffrestargs &optional value-symbol)
+(defmacro maref-function-picker
+    (type-symbol category ffrestargs &optional value-symbol)
   "Generate sexp to select on the various gsl_{vector,matrix}*_{get,set} functions."
   (cons 'cond
 	(mapcar (lambda (tp)
@@ -103,18 +103,18 @@
 		     ,@(if value-symbol
 			   (list (cl-cffi tp) value-symbol)
 			   (list (cl-cffi tp))))))
-		*array-element-types*))))
+		*array-element-types*)))
 
 (defmethod maref
     ((pointer #.+foreign-pointer-class+) index &optional index2
      (type 'double-float))
   (if index2
-      #.(maref-function-picker
-	 'type 'matrix
-	 '(:pointer pointer sizet index sizet index2))
-      #.(maref-function-picker
-	 'type 'vector
-	 '(:pointer pointer sizet index))))
+      (maref-function-picker
+	 type matrix
+	 (:pointer pointer sizet index sizet index2))
+      (maref-function-picker
+	 type vector
+	 (:pointer pointer sizet index))))
 
 ;;; Index the GSL function names to maref
 #.(cons 'progn
@@ -138,12 +138,12 @@
     (value (pointer #.+foreign-pointer-class+) index &optional index2
      (type 'double-float))
   (if index2
-      #.(maref-function-picker
-	 'type 'matrix
-	 '(:pointer pointer sizet index sizet index2) 'value)
-      #.(maref-function-picker
-	 'type 'vector
-	 '(:pointer pointer sizet index) 'value)))
+      (maref-function-picker
+       type matrix
+       (:pointer pointer sizet index sizet index2) value)
+      (maref-function-picker
+       type vector
+       (:pointer pointer sizet index) value)))
 
 ;;; Index the GSL function names to (setf maref)
 #.(cons 'progn
diff --git a/gsll-tests.asd b/gsll-tests.asd
index b025db793b75fb44154ceaa69ffdfcdc8aadc85f..8c7c0023751bd89d8f544625b83f6ed353a0f739 100644
--- a/gsll-tests.asd
+++ b/gsll-tests.asd
@@ -1,6 +1,6 @@
 ;; Definition of GSLL system 
 ;; Liam Healy
-;; Time-stamp: <2009-01-04 11:26:56EST gsll-tests.asd>
+;; Time-stamp: <2009-01-18 19:02:18EST gsll-tests.asd>
 ;; $Id$
 
 (asdf:defsystem "gsll-tests"
@@ -120,6 +120,7 @@
 	     (:file "negative-binomial")
 	     (:file "numerical-differentiation")
 	     (:file "numerical-integration")
+	     (:file "ode")
 	     (:file "pareto")
 	     (:file "permutation")
 	     (:file "poisson")
diff --git a/init/callback.lisp b/init/callback.lisp
index 8b52f56ab14cda1385b072c5ddfd39366a148140..32d333ba5b569565b25ca2a2716a5c59f9598a7e 100644
--- a/init/callback.lisp
+++ b/init/callback.lisp
@@ -1,6 +1,6 @@
 ;; Foreign callback functions.               
 ;; Liam Healy 
-;; Time-stamp: <2008-02-16 10:52:36EST callback.lisp>
+;; Time-stamp: <2009-01-18 17:04:46EST callback.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -61,17 +61,73 @@
   (parameters :pointer))
 
 ;;;;****************************************************************************
-;;;; Macros for defining a callback and placing in a structure
+;;;; Macros for defining a callback to wrap a CL function
 ;;;;****************************************************************************
 
 ;;; Usage example for scalar function (e.g. numerical-integration,
 ;;; numerical-differentiation, chebyshev, ntuple).  
 ;;; (defmcallback myfn :double :double)
-;;; Usage example for gsl-vector function (e.g. roots-multi)
+;;; Usage example for vector function (e.g. roots-multi)
 ;;; (defmcallback myfn :pointer :int (:pointer))
 ;;; Usage example for function and derivative
 ;;; (defmcallback fdf :pointer :double (:pointer :pointer))
 ;;; (defmcallback fdf :success-failure :int (:pointer :pointer))
+;;; Usage example for def-ode-functions
+;;; (defmcallback vanderpol :success-failure (:double (:double 2) (:set :double 2)))
+;;; or
+;;; (defmcallback vanderpol :success-failure (:double :pointer :pointer))
+;;; to read and set within the CL function with #'dcref.
+
+;;; (callback-args '(:double (:double 2) (:set :double 2)))
+;;; ((#:ARG1193 :DOUBLE) (#:ARG1194 :POINTER) (#:ARG1195 :POINTER))
+(defun callback-args (types)
+  "The arguments passed by GSL to the callback function."
+  (mapcar (lambda (type)
+	    (let ((symbol (gensym "ARG")))
+	      (list symbol
+		    (if (listp type)	; like (:double 3)
+			:pointer	; C array
+			type))))
+	  (if (listp types) types (list types))))
+
+;;; (embedded-clfunc-args '(:double (:double 2) (:set :double 2)) (callback-args '(:double (:double 2) (:set :double 2))))
+;;; (#:ARG1244 (MEM-AREF #:ARG1245 ':DOUBLE 0) (MEM-AREF #:ARG1245 ':DOUBLE 1))
+(defun embedded-clfunc-args (types callback-args)
+  "The arguments passed to the CL function call embedded in the callback."
+  (loop for spec in types
+     for (symbol type) in callback-args
+     append
+     (unless (and (listp spec) (eq (first spec) :set))
+       (if (listp spec)
+	   (loop for ind from 0 below (second spec)
+	      collect `(cffi:mem-aref ,symbol ',(first spec) ,ind))
+	   (list symbol)))))
+
+(defun callback-set-mvb (form types callback-args)
+  "Create the multiple-value-bind form in the callback to set the return C arrays."
+  (multiple-value-bind (settype setcba)
+      (loop for cba in callback-args
+	 for type in types
+	 for setting = (and (listp type) (eq (first type) :set))
+	 when setting
+	 collect cba into setcba
+	 when setting
+	 collect type into settype
+	 finally (return (values (mapcar 'rest settype) setcba)))
+    (print settype)
+    (let* ((setvbls (embedded-clfunc-args settype setcba))
+	   (count (apply '+ (mapcar 'second settype)))
+	   (mvbvbls (loop repeat count collect (gensym "SETCB"))))
+      (if (zerop count)
+	  form
+	  `(multiple-value-bind ,mvbvbls
+	       ,form
+	     (setf ,@(loop for mvbvbl in mvbvbls
+			for setvbl in setvbls
+			append (list setvbl mvbvbl))))))))
+
+;;; (DEFMCALLBACK VANDERPOL :SUCCESS-FAILURE (:DOUBLE (:DOUBLE 2) (:SET :DOUBLE 2)))
+;;; (DEFMCALLBACK VANDERPOL :SUCCESS-FAILURE (:DOUBLE (:DOUBLE 2)))
 
 (defmacro defmcallback
     (name &optional (return-type :double) (argument-types :double)
@@ -81,33 +137,39 @@
    argument-types is a single type or list of types of the argument(s)
    that appear before parameters, and the additional-argument-types
    (default none) is a single type or list of types of the argument(s)
-   that appear after parameters.  The return-type is the type that
-   should be returned to GSL.  If :success-failure, a GSL_SUCCESS
-   code (0) is always returned; if :pointer, a null pointer is
-   returned."
-  (flet ((arg-type (types)
-	   (when types
-	     (mapcar (lambda (type) (list (gensym "ARG") type))
-		     (if (listp types) types (list types))))))
-    (let ((arguments (arg-type argument-types))
-	  (additional-arguments (arg-type additional-argument-types)))
-      `(cffi:defcallback ,name
-	,(if (eq return-type :success-failure) :int return-type)
-	(,@arguments (params :pointer) ,@additional-arguments)
-	;; Parameters as C argument are always ignored, because we have
-	;; CL specials to do the same job.
-	(declare (ignore params))
-	(,name ,@(mapcar #'first (append arguments additional-arguments)))
-	,@(case
-	   return-type
-	   (:success-failure
-	    ;; We always return success, because if there was a
-	    ;; problem, a CL error would be signalled.
-	    '(success))
-	   (:pointer
-	    ;; For unclear reasons, some GSL functions want callbacks
-	    ;; to return a void pointer which is apparently meaningless.
-	    '((cffi:null-pointer))))))))
+   that appear after parameters.  The argument types are C types or
+   a list of a C type and a length, indicating a C array of that type
+   for which each element will be passed as a separate argument.
+   The return-type is the type that should be returned to GSL.
+   If :success-failure, a GSL_SUCCESS code (0) is always returned;
+   if :pointer, a null pointer is returned."
+  (let* ((atl (if (listp argument-types) argument-types (list argument-types)))
+	 (aatl (if (listp additional-argument-types) additional-argument-types
+		   (list additional-argument-types)))
+	 (cbargs (callback-args atl))
+	 (cbaddl (callback-args aatl)))
+    `(cffi:defcallback ,name
+	 ,(if (eq return-type :success-failure) :int return-type)
+	 (,@cbargs (params :pointer) ,@cbaddl)
+       ;; Parameters as C argument are always ignored, because we have
+       ;; CL specials to do the same job.
+       (declare (ignore params))
+       ,(callback-set-mvb
+	 `(,name
+	   ,@(append
+	      (embedded-clfunc-args atl cbargs) (embedded-clfunc-args aatl cbaddl)))
+	 (append atl aatl)
+	 (append cbargs cbaddl))
+       ,@(case
+	  return-type
+	  (:success-failure
+	   ;; We always return success, because if there was a
+	   ;; problem, a CL error would be signalled.
+	   '(success))
+	  (:pointer
+	   ;; For unclear reasons, some GSL functions want callbacks
+	   ;; to return a void pointer which is apparently meaningless.
+	   '((cffi:null-pointer)))))))
 
 (defmacro defcbstruct
     (functions &optional (structure 'gsl-function) additional-slots)
diff --git a/ordinary-differential-equations/ode-example.lisp b/ordinary-differential-equations/ode-example.lisp
index f347de2046998fcbb612634712322cb68e93eb6c..95a354b747f80f06d04fca308ad388e9e6b31bfc 100644
--- a/ordinary-differential-equations/ode-example.lisp
+++ b/ordinary-differential-equations/ode-example.lisp
@@ -1,45 +1,48 @@
 ;; Example ODE                               
 ;; Liam Healy Sat Sep 29 2007 - 17:49
-;; Time-stamp: <2008-01-28 22:21:07EST ode-example.lisp>
+;; Time-stamp: <2009-01-18 18:29:53EST ode-example.lisp>
 ;; $Id$
 
 ;;; van der Pol as given in Section 25.5 of the GSL manual.  To
 ;;; reproduce that example, (integrate-vanderpol 100.0d0)
 
+;;; Notice how the variable mu, which is a "parameter" in GSL's view,
+;;; to be passed on to the evaluation functions, is declared special
+;;; so its value is shared dynamically.  Thus GSL's parameters are
+;;; always empty (and not accessible through user functions in GSLL).
+
 (in-package :gsl)
 
-(defun vanderpol (time y dydt)
+(defun vanderpol (time y0 y1)
   (declare (special mu) (ignorable time))
-  (with-c-doubles ((y y0 y1) (dydt dydt0 dydt1))
-    (setf dydt0 y1
-	  dydt1 (- (- y0) (* mu y1 (- (* y0 y0) 1))))))
+  (values y1 (- (- y0) (* mu y1 (- (* y0 y0) 1)))))
 
-(defun vanderpol-jacobian (time y dfdy dfdt)
+(defun vanderpol-jacobian (time y0 y1)
   (declare (special mu) (ignorable time))
-  (with-c-doubles ((y y0 y1) (dfdy dfdy0 dfdy1 dfdy2 dfdy3) (dfdt dfdt0 dfdt1))
-    (setf dfdt0 0.0d0
-	  dfdt1 0.0d0
-	  dfdy0 0.0d0
-	  dfdy1 1.0d0
-	  dfdy2 (- (* -2 mu y0 y1) 1.0d0)
-	  dfdy3 (* -1 mu (- (* y0 y0) 1.0d0)))))
+  (values 0.0d0 0.0d0			; dfdt
+	  0.0d0 1.0d0			; dfdy in row-major order: 00, 01, 10, 11
+	  (- (* -2 mu y0 y1) 1.0d0)
+	  (* -1 mu (- (* y0 y0) 1.0d0))))
 
 (def-ode-functions vanderpol vanderpol-jacobian 2)
 
 (defparameter *max-iter* 2000)
 
-(defun integrate-vanderpol (max-time &optional (step-size 1.0d-6))
+(defun integrate-vanderpol (max-time &optional (step-size 1.0d-6) (print-steps t))
   "Integrate the van der Pol oscillator as given in Section 25.5 of the
    GSL manual.  To reproduce that example, (integrate-vanderpol 100.0d0)."
   (let ((mu 10.0d0) (time 0.0d0) (iter 0))
     (declare (special mu))
     (with-ode-integration (time step-size (dependent dep0 dep1) 2)
       (setf dep0 1.0d0 dep1 0.0d0)
-      (loop (when (or (>= (dcref time) max-time) (> iter *max-iter*)) (return))
-	    (apply-evolution
-	     evolve control stepper vanderpol
-	     time max-time step-size dependent)
-	    (incf iter)
-	    (format t "~&~12,6f~10t~12,6f~24t~12,6f"
-		    (dcref time) dep0 dep1)))))
-
+      (loop (when (or (>= (dcref time) max-time) (> iter *max-iter*))
+	      (return (values (dcref time) dep0 dep1)))
+	 (apply-evolution
+	  evolve control stepper vanderpol
+	  time max-time step-size dependent)
+	 (incf iter)
+	 (when print-steps
+	   (format t "~12,6f~10t~12,6f~24t~12,6f~&"
+		   (dcref time) dep0 dep1))))))
+
+(save-test ode (integrate-vanderpol 1.0d0 1.d-4 nil))
diff --git a/ordinary-differential-equations/ode-system.lisp b/ordinary-differential-equations/ode-system.lisp
index 04cceaea277762fe7d8f1cf54f42392f1da4c9d7..dafbbd3076c8d24fea51e21ae565847f9e4b7592 100644
--- a/ordinary-differential-equations/ode-system.lisp
+++ b/ordinary-differential-equations/ode-system.lisp
@@ -1,6 +1,6 @@
 ;; ODE system setup
 ;; Liam Healy, Sun Apr 15 2007 - 14:19
-;; Time-stamp: <2008-08-21 22:04:31EDT ode-system.lisp>
+;; Time-stamp: <2009-01-18 18:23:07EST ode-system.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -16,19 +16,29 @@
 (export '(def-ode-functions with-ode-integration))
 
 (defmacro def-ode-functions (name jacobian dimension)
-  "Setup functions for ODE integrators.
+  "Setup functions for ODE integrators.  The variable name is used as the name of the 
    The CL functions name and jacobian should be defined previously
-   with defuns."
-  ;; The function should take three arguments: time, dependent, derivatives
-  ;; The latter two will be C arrays.  To reference them, use #'with-c-doubles.
-  ;; To make this more transparent using a normal CL function
-  ;; would require transferring numbers back and forth between C and CL arrays,
-  ;; which could be inefficient.
+   with defuns.
+   The function and Jacobian arguments are the same: time and
+   dimension dependent variables as scalars.
+   The function returns dimension multiple values corresponding to
+   the derivatives of the dependent variables, f(t).
+   The Jacobian returns dimension^2 values corresponding to the
+   partial derivatives of each of the functions f with respect to each
+   of the dependent variables y, and dimension values corresponding to
+   the derivatives of f with respect to t."
+  ;; set return values
+  ;; Possible future improvements: take/set arrays easily, allow lambdas instead of named functions. 
   `(progn
-    (defmcallback ,name :success-failure (:double :pointer :pointer))
-    (defmcallback ,jacobian :success-failure (:double :pointer :pointer :pointer))
-    (defcbstruct (,name function ,jacobian jacobian) ode-system
-      ((dimension ,dimension)))))
+     (defmcallback
+	 ,name :success-failure
+       (:double (:double ,dimension) (:set :double ,dimension)))
+     (defmcallback
+	 ,jacobian :success-failure
+       (:double (:double ,dimension) (:set :double ,(expt dimension 2))
+		(:set :double ,dimension)))
+     (defcbstruct (,name function ,jacobian jacobian) ode-system
+       ((dimension ,dimension)))))
 
 (defmacro with-ode-integration
     ((time step-size dependent dimensions &optional (stepper '*step-rk8pd*)
diff --git a/tests/ode.lisp b/tests/ode.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..2a760d4f706db69afdd493b685e5d1bd591d1f87
--- /dev/null
+++ b/tests/ode.lisp
@@ -0,0 +1,11 @@
+;; Regression test ODE for GSLL, automatically generated
+
+(in-package :gsl)
+
+(LISP-UNIT:DEFINE-TEST ODE
+                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+                        (LIST 1.0d0 -1.4568622636249005d0
+                              -11.547385179410822d0)
+                        (MULTIPLE-VALUE-LIST
+                         (INTEGRATE-VANDERPOL 1.0d0 1.d-4 NIL))))
+