(defun Phi (phiE phiW source dx2 D L2)
  (/
    (+
      (/ (+ phiW phiE) dx2)
      (/ source D)
    ) 
    (+
      (/ 2 dx2)
      (/ 1 L2)
    )
  )
)

(defun PhiEdge (phiOpp albedo source dx2 D L2)
  (/
    (+
      (/
        (* 2 phiOpp)
        dx2
      )
      (/ source D)
    )
    (+
      (/ 2 dx2)
      (/
        (- 1 albedo)
        (* (+ 1 albedo) D (sqrt dx2))    
      )
      (/ 1 L2)
    )
  )
)

(defun FiniteDiff (flux source albedoR dx2 D L2)
  (cond
    ((null (cddr flux)) 
      (cons
        (car flux)
        (cons 
          (PhiEdge 
            (car flux)
            albedoR 
            (car source) 
            dx2 D L2
          )
          '()
        )
      )
    )
    (t
      (cons
        (car flux)
        (FiniteDiff
          (cons
            (Phi
              (caddr flux)
              (car flux)
              (car source)
              dx2 D L2
            )
            (cddr flux)
          )
          (cdr source)
          albedoR dx2 D L2
        )
      )
    )
  )
)

(defun SolveFlux (flux source albedoL albedoR dx2 D L2)
  (loop for i from 1 to 100 do
    (setq
      flux
      (cons
        (PhiEdge (cadr flux) albedoL (car source) dx2 D L2)
        (cdr flux)
      )
    )
    (setq
      flux 
      (FiniteDiff flux source albedoR dx2 D L2)
    )
  )
  (write flux)
)

(setq source '())
(setq flux '())
(loop for i from 0 to 20 do
  (setq source
    (append source 
      (cons (cos (/ (* i 1.5705) 20)) '())
    )
  )
  (setq flux (cons 1 flux))
)

(setq dx (/ 10 20))
(setq dx2 (* dx dx))
(setq D 1)
(setq L2 (/ D 0.1))

(SolveFlux flux source 0 0 dx2 D L2)
(SolveFlux flux source .2 .2 dx2 D L2)
(SolveFlux flux source .2 .9 dx2 D L2)