(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.
|#
    (let ((result (make-vector (fx+ n 1)))
	  (L2-inner-product-operator (->Linear-finite-element-operator fes))
	  (B (->Linear-finite-element-operator fes))
	  (delta-t (fl/ T (fixnum->flonum n)))
	  (inverse-delta-t (fl/ (fixnum->flonum n) T)))
      
      (define (predictor i)
	(if (fx= i 1)
	    (vector-ref result 0)
	    (Vector-subtract (Vector-scale 2.0 (vector-ref result (fx- i 1)))
			     (vector-ref result (fx- i 2)))))
      
      (define (solve A b . predictor)
	(if (null? predictor)
	    (conjugate-gradient A
				b 
				(lambda (A preconditioner z p r x c m)
				  (or (fl< c 1.0e-20)
				      (fx= m 100))))
	    (let* ((predictor (car predictor))
		   (residual (Vector-subtract b (Operator-apply A predictor)))
		   (increment (solve A residual)))
	      (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 (fx> backward-difference-steps 0)
	  (let ((backward-LHS-operator (Vector-add (Vector-scale inverse-delta-t L2-inner-product-operator) B)))
	    (do ((i 1 (fx+ i 1)))
		((fx> i (fxmin n backward-difference-steps)))
	      (let ((rhs-vector (instantiate Linear-finite-element-vector
					     space: fes))
		    (t^i (fl* (fixnum->flonum 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 (fx- i 1)))))
					     (predictor i)))))))
      (if (fx> 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 (fx+ backward-difference-steps 1) (fx+ i 1)))
		((fx> i n))
	      (let ((t^i-1/2 (fl* delta-t (fl- (fixnum->flonum 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 (fx- i 1)))
							 rhs-vector)
					     (predictor i)))))))
      result))