; Solve inhomogeneous TriDiagonal linear equations.
(define (solve aa ab ac b) (let* ((r (/ (car ab)))
      (p (* r (car b))) (q (- (* r (car ac)))))
   (let ((r1
     (let s ((p p) (q q) (aa aa) (ab (cdr ab)) (ac (cdr ac)) (b (cdr b)))
; (write (list "I am" p "+ you*" q)) (newline)
    (let ((r (/ (+ (* (car aa) q) (car ab))))
          (mec (- (car b) (* (car aa) p))))
    (if (null? ac) (list (* r mec))
        (let ((rp (s (* mec r) (- (* (car ac) r))
             (cdr aa) (cdr ab) (cdr ac) (cdr b))))
;        (write (list "we are" rp)) (newline)
         (cons (* r (- mec (* (car ac) (car rp)))) rp)))))))
   (cons (+ p (* q (car r1))) r1))))

(define (ml aa ab ac x)
  (cons (+ (* (car ab) (car x)) (* (car ac) (cadr x)))
  (let p ((aa aa)(ab (cdr ab)) (ac (cdr ac)) (x x))
   (let ((tw (+ (* (car aa) (car x)) (* (car ab) (cadr x)))))
    (if (null? ac) (list tw)
       (cons (+ tw (* (car ac) (caddr x)))
          (p (cdr aa) (cdr ab) (cdr ac) (cdr x))))))))

(define (pf c) (apply (lambda (aa ab ac x)
(write (list "aa ab ac x is " aa ab ac x)) (newline)
(let* ((b (ml aa ab ac x))
      (xz (solve aa ab ac b)))
(write (list "b is " b)) (newline)
(write (list "x is " x)) (newline)
(write (list "xz is" xz)) (newline)
    (let r ((x x) (xz xz)) (if (null? x) 0
      (max (abs (- (car x) (car xz))) (r (cdr x) (cdr xz)))))))
c))

(define (I) (pf '((0 0 0)
     (1 1 1 1)
     (0 0 0)
     (2 3 5 7))))

(define (I) (pf '((0.21 0.1 0.16)
     (1.4 0.3 1.6 -1.44)
     (0.21 0.18 0.14)
     (2 6 5 7))))

ab0*x0 + ac0*x1 = b0
x0 = (b0 - ac0*x1)/ab0 = b0/ab0 + (-ac0/ab0)*x1 = p + q*x1

aa0*x0 + ab1*x1 + ac1*x2 = b1
aa0*(p + q*x1) + ab1*x1 + ac1*x2 = b1
aa0*p + (aa0*q + ab1)*x1 + ac1*x2 = b1
(aa0*q + ab1)*x1 + ac1*x2 = b1 - aa0*p
(aa0*q + ab1)*x1 = b1 + (-ac1)*x2
x1 = (b1 - aa0*p)/(aa0*q + ab1) + ((-ac1)/(aa0*q + ab1))*x2

aan*xn' + abn*xn = bn
aan*(p + q*xn) + abn*xn = bn
aan*p + aan*q*xn + abn*xn = bn
xn*(aan*q + abn) = bn - aan*p
xn = (bn - aan*p)/(aan*q + abn)