(require (planet "ode-initval.ss" ("williams" "science.plt" 2 0)))
(require (lib "plot.ss" "plot"))

(define (func t y f params)
  (let ((y0 (vector-ref y 0)))
    (vector-set!
     f 0
     (* (- 1500.0 (vector-ref y 0)) 0.12))))

(define (main)
  (let* ((type rk4-ode-type)
         (step (make-ode-step type 1))
         (control #f) ;(control-y-new 1.0e-6 0.0))
         (evolve (make-ode-evolve 1))
         (system (make-ode-system func #f 1 '()))
         (t (box 0.0))
         (t1 20.0)
         (h (box (/ 1.0 60.0)))
         (y #(150.0))
         (y0-values '()))
    (let loop ()
      (if (< (unbox t) t1)
          (begin
            (ode-evolve-apply
             evolve control step system
             t t1 h y)
            ;(printf "~a ~a ~a~n" (unbox t) (vector-ref y 0) (vector-ref y 1))
            (set! y0-values (cons (vector (unbox t) (vector-ref y 0)) y0-values))
            (loop))))
    (printf "Number of iterations   = ~a~n" (ode-evolve-count evolve))
    (printf "Number of failed steps = ~a~n" (ode-evolve-failed-steps evolve))
    (printf "~a~n" (plot (points (reverse y0-values))
                         (x-min 0.0)
                         (x-max 20.0)
                         (y-min 0.0)
                         (y-max 1500.0)))
    ))

(main)