We consider a fully practical finite element approximation of the nonlinear degenerate parabolic system gammapartial derivativeu/partial derivativet - del.(b(mu) del[w+alphaphi]) = 0, w = -gammaDeltau+gamma(-1)Psi'(u), del.(c(u) delphi) =0 subject to an initial condition u(0)(.) is an element of[-1, 1] on u and flux boundary conditions on all three equations. Here gamma is an element ofR(>0), alpha is an element ofR(greater than or equal to0), Psi is a nonsmooth double well potential, and c(u) := 1+u, b(u) := 1-u(2) are degenerate coefficients. The degeneracy in b restricts u(.,.) is an element of[-1, 1]. The above, in the limit gamma --> 0, models the evolution of voids by surface diffusion in an electrically conducting solid. In addition to showing stability bounds for our approximation, we prove convergence, and hence existence of a solution to this nonlinear degenerate parabolic system in two space dimensions. Furthermore, an iterative scheme for solving the resulting nonlinear discrete system is introduced and analyzed. Finally, some numerical experiments are presented.