; Apply an iterator, such as RK4 or Euler

(defun apply-iterator (iterator f xo to tf dt)
  (do ((time to (+ time dt))
       (x xo (funcall iterator f x time dt)))
      ((> time tf) (format t "~%"))
    (format t "~A       ~A~%" time x)))

; Runge-Kutta 4th order iterator

(defun rk4 (f xo to dt)
    (let* ((k1 (funcall f xo to))
           (k2 (funcall f (+ xo (* k1 dt 0.5)) (+ to (/ dt 2))))
           (k3 (funcall f (+ xo (* k2 dt 0.5)) (+ to (/ dt 2))))
           (k4 (funcall f (+ xo (* k3 dt)) (+ to dt))))
      (+ xo (* (/ dt 6) (+ k1 (* 2 k2) (* 2 k3) k4)))))

; Euler iterator

(defun euler (f xo to dt)
  (+ xo (* dt (funcall f xo to))))

(format t "First exercise: p246, example 3~%~%")

(apply-iterator #'rk4 
                #'(lambda (xo to) (+ to 1 (- xo))) 
                1.0 0.0 1.01 0.1)

(format t "Second exercise: p249, problem 7a (euler and rk4)~%~%")

(defun f7a (x time)
  (- (expt time -2) (/ x time)))

(format t "apply Euler~%")
(apply-iterator #'euler #'f7a -1.0 1.0 2.01 0.025)

(format t "apply RK4~%")
(apply-iterator #'rk4 #'f7a -1.0 1.0 2.01 0.1) 

(format t "correct answer (which is better?)~%")
(apply-iterator #'rk4 #'f7a -1.0 1.0 2.01 0.01)