BCa = [0;0] ;
BCb = [1;2] ;

N = [ 10, 20, 40, 80, 160, 320 ] ;

for i=1:length(N)
  [x,y] = BVPsolver( @p1, @q1, @r1, BCa, BCb, N(i) );
  err(i) = norm( y - exact1( x ), inf ) ;
end

for i=2:length(N)
  h1 = 1/N(i-1) ;
  h  = 1/N(i) ;
  p  = log( err(i-1)/err(i) ) / log( h1/h ) ;
  disp(p) ;
end
