N = 5 ;
[ A, b, X1, Y1, Z1 ] = elliptic( N, @f, @g ) ;
N = N*2
[ A, b, X2, Y2, Z2 ] = elliptic( N, @f, @g ) ;
N = N*2
[ A, b, X3, Y3, Z3 ] = elliptic( N, @f, @g ) ;
N = N*2
[ A, b, X4, Y4, Z4 ] = elliptic( N, @f, @g ) ;
N = N*2
[ A, b, X5, Y5, Z5 ] = elliptic( N, @f, @g ) ;

EX1 = g(X1, Y1) ;
EX2 = g(X2, Y2) ;
EX3 = g(X3, Y3) ;
EX4 = g(X4, Y4) ;
EX5 = g(X5, Y5) ;

E1 = max(max( abs(Z1 - EX1) )) ;
E2 = max(max( abs(Z2 - EX2) )) ;
E3 = max(max( abs(Z3 - EX3) )) ;
E4 = max(max( abs(Z4 - EX4) )) ;
E5 = max(max( abs(Z5 - EX5) )) ;

disp( sprintf( 'E1 = %g', E1 ) ) ; 
disp( sprintf( 'E2 = %g', E2 ) ) ; 
disp( sprintf( 'E3 = %g', E3 ) ) ; 
disp( sprintf( 'E4 = %g', E4 ) ) ; 
disp( sprintf( 'E5 = %g', E5 ) ) ; 

disp( sprintf( 'order = %g', log( E1/E2 ) / log( 2 ) ) ) ;
disp( sprintf( 'order = %g', log( E2/E3 ) / log( 2 ) ) ) ;
disp( sprintf( 'order = %g', log( E3/E4 ) / log( 2 ) ) ) ;
disp( sprintf( 'order = %g', log( E4/E5 ) / log( 2 ) ) ) ;
