(* This is the second-order symbolic solver we worked through, and its check. *)

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 } ;

v = LinearSolve [ a , c ] ;
p = p1 * e ; q = q1 * e ; r = r1 * e ;
rule1 = e ^ k_ /; k > 2 -> 0 ;
w = ExpandAll [ Together [ v ] ] /. rule1 ;

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 )
] ;

rule2 = k1_ / k2_ :> k1 * invert [ k2 ] ;
x = Expand [ w //. rule2 ] /. rule1 ;
p =. ; q =. ; r =. ;
p1 = p / e ; q1 = q / e ; r1 = r / e ;
y = x ;
p1 =. ; q1 =. ; r1 =. ;
Print [ y ] ;

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 ] ;
t /. x_ /; Abs[x] < 1.0*^-11 -> 0