(* This is the Gaussian Elimination symbolic solver, and its check, made
second-order; see the file gaussian_4.math for comments on its basic
construction. *)

a = { {  4.2 ,      2.2 + p , -3.9 ,      9.3 ,      0.1     } ,
      {  8.6 - p ,  r ,        0.7 ,     -2.3 ,     -0.3     } ,
      {  8.4 ,     -5.9 ,     -8.1 + q ,  9.6 ,      3.8     } ,
      { -0.8 ,     -9.4 ,     -9.9 ,      9.9 ,      5.0 - r } ,
      { -1.3 ,     -8.1 ,      0.6 ,     -9.2 + r , -7.3     } } ;

c = { p , 1.0 , p + q , 1.0 , q } ;

(* We use Expand and apply rule1 to clobber high-order terms following
every multiplication - they obviously can't be created by an addition or
subtraction. *)

rule1 = e ^ k_ /; k > 2 -> 0 ;

(* We call invert to change each division into a multiplication, to
second-order.  It could easily be changed to be accurate to Nth order,
for any suitable N. *)

invert [ arg_ ] := Module [ { x, y , z , n } ,
    x = CoefficientList [ arg , e ] ;
    n = Dimensions [ x ] [[ 1 ]] ;
    If [ n < 2 , ( y = 0.0 ; z = 0.0 ) ,
        If [ n < 3 ,( y = x [[ 2 ]] / x [[ 1 ]] ; z = 0.0 ) ,
            ( y = x [[ 2 ]] / x [[ 1 ]] ; z = x [[ 3 ]] / x [[ 1 ]] )
        ]
    ] ;
    (1.0 / x [[ 1 ]] ) * ( 1.0 - y * e + ( y ^ 2 - z ) * e ^ 2 )
] ;

(* The function symsolver is almost unchanged, except that each
division and application of t2 is replaced by a call to invert. *)

symsolver [ a_ , c_ , t1_ ] := Module [ { a1 , c1 , i, j, n, x } ,
    a1 = a ; c1 = c ;
    n = Dimensions [ a1 ] [[ 1 ]] ;
    p = p1 * e ; q = q1 * e ; r = t1 * e ;
    For [ j = 1 , j < n , ++ j ,
        For [ i = j + 1 , i <= n , ++ i ,
            x = Expand [ a1 [[ i , j ]] *
                invert [ a1 [[ j , j ]] ] ] /. t1 ;
            a1 [[ i , ;; ]] = a1 [[ i , ;; ]] -
                Expand [ x * a1 [[ j , ;; ]] ] /. t1 ;
            c1 [[ i ]] = c1 [[ i ]] -
                Expand [ x * c1 [[ j ]] ] /. t1 ;
        ]
    ] ;
    For [ j = n , j > 1 , -- j ,
        For [ i = 1 , i < j , ++ i ,
            x = Expand [ a1 [[ i , j ]] *
                invert [ a1 [[ j , j ]] ] ] /. t1 ;
            a1 [[ i , ;; ]] = a1 [[ i , ;; ]] -
                Expand [ x * a1 [[ j , ;; ]] ] /. t1 ;
            c1 [[ i ]] = c1 [[ i ]] -
                Expand [ x * c1 [[ j ]] ] /. t1 ;
        ]
    ] ;
    For [ i = 1 , i <= n , ++ i ,
        c1 [[ i ]] *= invert [ Diagonal [ a1 ] [[ i ]] ]
    ] ;
    (* The need to clobber multiples of zero explicitly is yet
    another artifact of Mathematica's processing - sigh. *)
    Expand [ c1 ] /. t1 /. 0.0 * _ -> 0.0
] ;

(* We convert p, q and r to multiples of the infinitesimal, e, apply the
solver, reverse the conversions, and print the result converted back to
the original variables. *)

p = p1 * e ; q = q1 * e ; r = r1 * e ;
x = symsolver [ a , c , rule1 ] ;
p =. ; q =. ; r =. ;
p1 = p / e ; q1 = q / e ; r1 = r / e ;
y = x ;
p1 =. ; q1 =. ; r1 =. ;
Print [ y ] ;

(* We do exactly the same process, but print the solution multiplied by
the matrix to check that we get back to the original vector, rather than
solving the equations.  We need to use Expand and rule1 again, as we may
reintroduce some second-order terms by doing this. *)

p = p1 * e ; q = q1 * e ; r = r1 * e ;
s = Expand [ a . y ] /. rule1 ;
p =. ; q =. ; r =. ;
p1 = p / e ; q1 = q / e ; r1 = r / e ;
t = s ;
p1 =. ; q1 =. ; r1 =. ;
Print [ t ] ;

(* We print the result with tiny floating-point values converted to
zero, because it clarifies the result; this is not something that should
be done if you are going to reuse the result. *)

t /. x_ /; Abs[x] < 1.0*^-11 -> 0