At steady state, the dissociation-association electrochemical reactions in a dilute solution are governed by the nondimensional equations
\[
∂
∂
x
(
d
i
∂
u
i
∂
x
+
d
i
e
i
∂
ϕ
∂
x
u
i
)
+
R
i
=
0
,
i
=
1
,
.
.
.
,
m
,
0
≤
x
≤
1
,
∑
i
=
1
m
e
i
u
i
=
0
\frac {\partial }{{\partial x}}\left ( {{d_i}\frac {{\partial {u_i}}}{{\partial x}} + {d_i}{e_i}\frac {{\partial \phi }}{{\partial x}}{u_i}} \right ) + {R_i} = 0, \qquad i = 1,...,m, 0 \le x \le 1, \\ \sum \limits _{i = 1}^m {{e_i}{u_i} = 0}
\]
where
d
i
,
e
i
{d_i}, {e_i}
, and
u
i
{u_i}
are the diffusion coefficient, charge, and concentration of the
i
i
th species, respectively.
ϕ
\phi
is the electric potential of the solution.
R
i
{R_i}
is the net source of the
i
i
th species in the solution which includes the external input or production due to the reactions in the solution. The extra electroneutrality condition
∑
i
=
1
m
e
i
u
i
=
0
\sum _{i = 1}^m {e_i}{u_i} = 0
determines the electric potential
ϕ
\phi
. This system of nonlinear differential equations is subject to the nonlinear boundary conditions modeling the actual electrode kinetics. We prove the existence of the solution to this system and construct an iterative numerical algorithm to compute the solution. Numerical results are also presented in the paper.