Cálculo de la curva de titulacion de un ácido poliprótico

\[ V_{b}\left[\frac {C_{b} K_{bc} [H^+] }{K_{bc}[H^+]+ K_{we} } +[H^+]- \frac{K_{we}} {[H^+]} \right ] = V_a \left[ C_a \left( \frac{K_{1c}[H^+]^2+2K_{1c}K_{2c}[H^+]+3K_{1c}K_{2c}K_{3c}}{[H^+]^3+ K_{1c}[H^+]^2+K_{1c}K_{2c}[H^+]+K_{1c}K_{2c}K_{3c}} \right )- [H^+] + \frac{K_{we}}{[H^+]} \right ] \]

## De acido 
Va=20 

## Concentraciones de la base y del acido
Cb= 0.1
Ca= 0.1

## Pka 
pKa1= 2.12
pKa2= 7.21
pKa3= 12.67

## Constantes de disociación de ácido
K1= 10^(-pKa1)
K2= 10^(-pKa2)
K3= 10^(-pKa3)

## Constantes de disociación del agua
Kwe= 1*10^-14

## Constante de disosiación basica
Kbc= 1 

pH=seq(1,12)
  
H= 10^-pH
Vb= ((Va)*( Ca * ( ( (K1*H^2) + (2*K1*K2*H) + (3*K1*K2*K3) ) / ( (H^3)+ (K1*H^2) + (K1*K2*H) + (K1*K2*K3) ) ) - (H) + (Kwe/H) )) / ( ((Cb*Kbc*H)/( (Kbc*H) + Kwe )) + (H) - (Kwe/H) )
print(Vb)

##  [1] -9.294908  6.024778 17.298794 19.712389 20.092225 21.158534 27.628121
##  [8] 37.210039 39.691454 40.074366 41.067097 51.143136