  In this paper we consider a coupled bulk-surface PDE in two space dimensions.
The model consists of a PDE in the bulk that is coupled to another PDE on the
surface through general nonlinear boundary conditions. For such a system we
propose a novel method, based on coupling a virtual element method [Beir\~ao da
Veiga et al., 2013] in the bulk domain to a surface finite element method
[Dziuk & Elliott, 2013] on the surface. The proposed method, which we coin the
Bulk-Surface Virtual Element Method (BSVEM) includes, as a special case, the
bulk-surface finite element method (BSFEM) on triangular meshes [Madzvamuse &
Chung, 2016]. The method exhibits second-order convergence in space, provided
the exact solution is $H^{2+1/4}$ in the bulk and $H^2$ on the surface, where
the additional $\frac{1}{4}$ is required only in the simultaneous presence of
surface curvature and non-triangular elements. Two novel techniques introduced
in our analysis are (i) an $L^2$-preserving inverse trace operator for the
analysis of boundary conditions and (ii) the Sobolev extension as a replacement
of the lifting operator [Elliott & Ranner, 2013] for sufficiently smooth exact
solutions. The generality of the polygonal mesh can be exploited to optimize
the computational time of matrix assembly. The method takes an optimised
matrix-vector form that also simplifies the known special case of BSFEM on
triangular meshes [Madzvamuse & Chung, 2016]. Three numerical examples
illustrate our findings.
