Michael Whiston (2071, Spring 13) asked a question about convergence
in Lab 4, Exercise 9, where the solution to the Dirichlet BVP
y''+y'+y=x*(1-x) is found using quadratic Lagrange elements in 1D on a
uniform mesh.  Comparing with the exact solution gives the table
  N    h        error      ratio
  7  1.2500e-1 7.2202e-06 15.460
 15  6.2500e-2 4.6703e-07 15.788
 31  3.1250e-2 2.9581e-08 14.013
 61  1.6129e-2 2.1110e-09 14.951
121  8.1967e-3 1.4119e-10
Why is the progression of ratios not monotone?

NOTE: error is measured in the inf norm.  Measuring in the 2-norm,
which is proportional to the L^2 norm, gives a similar progression:
   2.2508e-05 15.8535   
   1.4198e-06 15.9629   
   8.8941e-08 14.0836   
   6.3152e-09 15.0180
   4.2051e-10

Part of this question is roundoff error.  The matrix has condition
number around 1e4 or 1e5, so N=121 is about as fine a mesh as you 
can get in double precision.  I have looked at finding the solution
when the exact solution is x^3 and found an equally strange-looking
progression.  And why, after all, is this O(h^4)?  Shouldn't it be
O(h^3)?  Is that a function of the uniform mesh?

My plan is to re-write the Matlab as Fortran and use quad precision.
I should be able to get a few more mesh doublings that way.  I will
also handle the mesh so that it can be non-uniform and put that
question to bed.

1. Get the Matlab/octave version running, then start converting code.
* tested with testit.f90, testsolver.f90, and case soln=x*(1-x) in
  solve9
* ran solve9:
  Results: relative 2-norm
   N=           7  err=  2.25082372698789977944996516368153355E-0005
   N=          15  err=  1.41976120698713675663303820471654439E-0006
   N=          31  err=  8.89411404844063349763423227504439248E-0008
   N=          61  err=  6.31513613271189600494497534019060715E-0009
   N=         121  err=  4.21285661993855810578875401384659154E-0010
   N=         241  err=  2.72127114432575728708283198904980239E-0011
   N=         481  err=  1.72921713890222200081758987213504322E-0012
   N=         961  err=  1.08977918316003134947878547524348258E-0013
   N=        1921  err=  6.83951870934140758104742070770058018E-0015
  539.493u 0.680s 9:01.96 99.6%   0+0k 0+16io 0pf+0w
   N=        3841  err=  4.28360782372179875910273306266498229E-0016
  5673.030u 2.840s 1:34:54.79 99.6%       0+0k 0+16io 0pf+0w
* resulting ratios are
     15.854
     15.963
     14.084
     14.990
     15.481
     15.737
     15.868
     15.934
     15.967

2. Change solver to handle banded matrices.  Retain full matrix storage.
Visually compared errors against above: agree.  Time for 3841 case=2.67 min
on oz125
   N=           7  err=  2.25082372698789977944996516368153355E-0005
   N=          15  err=  1.41976120698713675663303820471654439E-0006
   N=        1921  err=  6.83951870934140758104742070770058018E-0015
   N=        3841  err=  4.28360782372179875910273306266498229E-0016
   N=        7681  err=  2.68004418722241950969193500373076638E-0017
   ratio 15.983