(* This is the 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 > 1 -> 0 ;
w = ExpandAll [ Together [ v ] ] /. rule1 ;
rule2 = k1_ / ( k2_ + k3_ * e ) ->
    k1 * ( 1 - ( k3 / k2 ) * e ) / k2 ;
x = w //. rule2 ;
y = Expand [ x ] /. rule1 ;
p =. ; q =. ; r =. ;
z = y /. e p1 -> p /. e q1 -> q /. e r1 -> r ;
Print [ z ] ;

p = p1 * e ; q = q1 * e ; r = r1 * e ;
s = Expand [ a . z ] /. rule1 ;
p =. ; q =. ; r =. ;
t = s /. e p1 -> p /. e q1 -> q /. e r1 -> r
Print [ t ] ;
t /. x_ /; Abs[x] < 1.0*^-12 -> 0