\[ 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(6)
  
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
n=3
pKa=c(2.1,7.21,12.67)
ka=10**(-pKa)

PH=seq(1,13,0.1)
H=10**(-PH)

D=H*n+ka[1]*H**(n-1)+ka[1]*ka[2]*H*(n-2)+ka[1]*ka[2]*ka[3]


alfa0=H**n/D
alfa1=ka[1]*H*(n-1)/D
alfa2=ka[1]*ka[2]*H*(n-2)/D
alfa3=ka[1]*ka[2]*ka[3]/D
plot(PH,alfa0,type="line")
## Warning in plot.xy(xy, type, ...): plot type 'line' will be truncated to first
## character