(define (parabolic-equation-solution a c f g gamma u0 fes T n backward-difference-steps)
#|
Solves the equation
$$
\alignedat 2
u_t-\nabla\cdot (a\nabla u) + cu&=f,&&\quad x\in\Omega,\ 0<t\leq T,
\\
a\frac{\partial u}{\partial n}+\gamma u&=g,&&\quad x\in\partial\Omega,\ 0<t\leq T,
\\
u(x,0)&=u_0(x),&&\quad x\in\Omega,
\endalignedat
$$
at times $t^i=i\Delta t$, $\Delta t=T/n$, $i=0,\dots,n$, in the linear finite
element space fes,
using backward-difference-steps backward-difference steps at the beginning
of the simulation and Crank-Nicholson steps for the rest of the simulation.

The linear solver is conjugate gradient without a preconditioner with a
ridiculously small tolerance; constant prediction is used for the
solution at time $t^1$, and linear prediction thereafter.
|#
  (declare (fixnum))
   
    (let ((result (make-vector (+ n 1)))
	  (L2-inner-product-operator (->Linear-finite-element-operator fes))
	  (B (->Linear-finite-element-operator fes))
	  (delta-t (FLOAT (/ T (##flonum.<-fixnum n))))
	  (inverse-delta-t (FLOAT (/ (##flonum.<-fixnum n) T))))
      
      (define (predictor i)
	(if (= i 1)
	    (vector-ref result 0)
	    (Vector-subtract (Vector-scale 2.0 (vector-ref result (- i 1)))
			     (vector-ref result (- i 2)))))
      
      (define (solve A b . predictor)
	(if (null? predictor)
	    (conjugate-gradient A
				b 
				(lambda (A preconditioner z p r x c m)
				  (or (FLOAT (< c 1.0e-20))
				      (FIX (= m 100)))))
	    (let* ((predictor (car predictor))
		   (residual (Vector-subtract b (Operator-apply A predictor)))
		   (increment (conjugate-gradient A
						  residual 
						  (lambda (A preconditioner z p r x c m)
						    (or (FLOAT (< c 1.0e-20))
							(FIX (= m 100)))))))
	      (Vector-add predictor increment))))
      
      (add-mass-terms! L2-inner-product-operator (lambda (p) 1.0))
      (add-stiffness-terms! B a)
      (add-mass-terms! B c)
      (add-boundary-terms! B gamma)
      
      ;;calculate u(0)
      (let ((rhs-vector (instantiate Linear-finite-element-vector
				     space: fes)))
	(add-interior-integral-terms! rhs-vector u0)
	(vector-set! result 0 (solve L2-inner-product-operator rhs-vector)))
      
      (if (> backward-difference-steps 0)
	  (let ((backward-LHS-operator (Vector-add (Vector-scale inverse-delta-t L2-inner-product-operator) B)))
	    (do ((i 1 (+ i 1)))
		((> i (min n backward-difference-steps)))
	      (let ((rhs-vector (instantiate Linear-finite-element-vector
					     space: fes))
		    (t^i (FLOAT (* (##flonum.<-fixnum i) delta-t))))
		(add-interior-integral-terms!  rhs-vector (lambda (p) (f t^i p)))
		(add-boundary-integral-terms!  rhs-vector (lambda (p) (g t^i p)))
		(vector-set! result i (solve backward-LHS-operator
					     (Vector-add rhs-vector 
							 (Operator-apply L2-inner-product-operator
									 (Vector-scale inverse-delta-t (vector-ref result (- i 1)))))
					     (predictor i)))))))
      (if (> n backward-difference-steps)
	  (let ((CN-LHS-operator (Vector-add (Vector-scale inverse-delta-t L2-inner-product-operator)
					     (Vector-scale 0.5 B)))
		(CN-RHS-operator (Vector-subtract (Vector-scale inverse-delta-t L2-inner-product-operator)
						  (Vector-scale 0.5 B))))
	    (do ((i (+ backward-difference-steps 1) (+ i 1)))
		((> i n))
	      (let ((t^i-1/2 (FLOAT (* delta-t (- (##flonum.<-fixnum i) 0.5))))
		    (rhs-vector (instantiate Linear-finite-element-vector
					     space: fes)))
		(add-interior-integral-terms!  rhs-vector (lambda (p) (f t^i-1/2 p)))
		(add-boundary-integral-terms!  rhs-vector (lambda (p) (g t^i-1/2 p)))
		(vector-set! result i (solve CN-LHS-operator
					     (Vector-add (Operator-apply CN-RHS-operator
									 (vector-ref result (- i 1)))
							 rhs-vector)
					     (predictor i)))))))
      result))