Hi Liam (& others) I am getting funny behavior with the ODE stepper: (As refresher, I am solving a 2nd order system that should reproduce a sinusoid) A trace of the dydt function (sin-ode below) shows that it is always called with arguments (time 0 0), even though the initial conditions are (time 0 1) (see the let form). It looks as if the initial conditions (stored in vector y below) are not passed as arguments to the dydt function (sin-ode). Or am I missing something? I tried both the rk2 and rk4 steppers and different step sizes. Thanks, Mirko Here is the source: ;;;; Example of a simple time-stepper (defun sin-ode (time y z) "Define ODE for a sinusoid y''=-y or as a system: y'=z z'=-y with initial conditions: y0=0 z0=1 " (declare (ignore time)) (values z (- y))) (defun sin-ode-jacobian (time y z) (declare (ignore time y z)) (values 0d0 0d0 0d0 1d0 -1d0 0d0)) #| (let ((time 0d0) (delta-t 0.01d0) (stepper-type +step-rk2+)) (let ((stepper (make-ode-stepper stepper-type 2 #'sin-ode #'sin-ode-jacobian)) (y (make-marray 'double-float :dimensions 2)) (dydt-in (make-marray 'double-float :dimensions 2)) (dydt-out (make-marray 'double-float :dimensions 2)) (yerr (make-marray 'double-float :dimensions 2))) (setf (maref y 0) 0d0 ;; y(0)=0 (maref y 1) 1d0);; y'(0)=1 (multiple-value-bind (yp zp) (sin-ode time (maref y 0) (maref y 1)) (setf (maref dydt-out 0) yp (maref dydt-out 1) zp)) (dotimes (i 2) (incf time delta-t) (setf (maref dydt-in 0) (maref dydt-out 0) (maref dydt-in 1) (maref dydt-out 1)) (format t "Y-in: ~12,4f, ~a, ~a~%" time (maref y 0) (maref y 1)) (format t "dydt-in:~12,4f, ~a, ~a~%" time (maref dydt-in 0) (maref dydt-in 1)) (apply-step stepper time y delta-t yerr dydt-in dydt-out) (format t "Y-out: ~12,4f, ~a, ~a~%" time (maref y 0) (maref y 1)) (format t "dydt-out:~12,4f, ~a, ~a~%~%" time (maref dydt-out 0) (maref dydt-out 1))) y)) |#