  We derive a numerical method, based on operator splitting, to abstract
parabolic semilinear boundary coupled systems. The method decouples the linear
components which describe the coupling and the dynamics in the bulk and on the
surface, and treats the nonlinear terms by approximating the integral in the
variation of constants formula. The convergence proof is based on estimates for
a recursive formulation of the error, using the parabolic smoothing property of
analytic semigroups and a careful comparison of the exact and approximate
flows. Numerical experiments, including problems with dynamic boundary
conditions, reporting on convergence rates are presented.
