mixture of multiple acids or multiple bases
principle : superposition ( = charge balance of species)
--------------------------------------------------------
--------------------------------------------------------
info:quantitatve  chemical analysis (fourth edition)
     I.M.KOLTHOFF,The Macmillan Company.
--------------------------------------------------------
peter.vlasschaert@gmail.com,30/09/2016

example :combination of 'two weak acids,both mono protic' 
         k(number component)A(number anion solution)
        example : k(1)A(1) = k1A1
first  : HA1=H(+) + A1(-) => k1a1 = ([H(+)]*[A1(-)])/[HA1] (1)
ca1    = [HA1]+[A1(-)] = ha1 +a1                           (2)
        example : k(2)A(2) = k2A2
second : HA2=H(+) + A1(-) => k2a2 = ([H(+)]*[A2(-)])/[HA2] (3)
ca2    = [HA2]+[A2(-)] = ha2 +a2                           (4)
species solution : negative  oh=[OH(-)],a1=[A1(-)],a2=[A2(-)]
                   positive  h=[H(+)]
charge balance   : [H(+)] = [OH(-)]+[A1(-)]+[A2(-)]
total  balance   : [H(+)]tot = [H(+)]'aq'+[H(+)]'HA1'+[H(+)]'HA2'
 result:           [H(+)]    = kw/[H(+)]+[A1(-)]+[A2(-)]
 math                h       = kw/h + a1 + a2              (5)

multiple bases:mono hydric bases : M1OH,M2OH,k1b1,k2b2,cb1,cb2
=>                                'mono hydric = mono hydroxy'
--------------------------------------------------------------
[OH(-)]tot = [OH(-)]'aq'+[OH(-)]'M1OH'+[OH(-)]'M2OH'

kill(all);

from equation (1) and (3)

p1:k1a1*ha1 = h*a1$
p2:k2a2*ha2 = h*a2$
p3:part(solve(p1,a1),1);
p4:part(solve(p2,a2),1);

from equation (5),put in a1,a2

p5:h=kw/h+rhs(p3)+rhs(p4)$

p6: lhs(p5)-rhs(p5)=0;

p7:ratsimp(lhs(p6))=0;

p8:num(lhs(p7))=0;

? only positive root has physical meaning
  'concentration = positive value.'

p9:part(solve(p8,h),2);

p10:ca1=ha1+a1;
p11:ca1=ha1+rhs(p3)$
p12:ca1=factor(rhs(p11));
p13:ca2=ha2+a2;
p14:ca2=ha2+rhs(p4)$
p15:ca2=factor(rhs(p14));

p16:part(solve(p12,ha1),1);
p17:part(solve(p15,ha2),1);

p18:ev(p9,p16,p17);

first approximation before '(p18)' => kw=0

p19:ev(p18,kw=0);

p20:ca1=0.01;
p21:ca2=0.1;
p22:k1a1=0.0001;
p23:k2a2=0.00001;

p24:ev(p19,p20,p21,p22,p23);

? h : for right hand side (rhs)
      => 'h>k1a1 and h>k2a2,kw=0'

? First approximation,h

p25:ha11=ev(rhs(p16),k1a1=0);
p26:ha22=ev(rhs(p17),k2a2=0);

p27:ev(p9,kw=0,ha1=rhs(p25),ha2=rhs(p26));

? second approximation,h

p28:h=ev(rhs(p24),p27);

p29:h=ev(rhs(p28),p20,p21,p22,p23);