We propose a new numerical method to solve the Cahn-Hilliard equation coupled with nonlinear wetting boundary conditions. We show that the method is mass-conservative and that the discrete solution satisfies a discrete energy law similar to the one satisfied by the exact solution. We perform several tests inspired by realistic situations to verify the accuracy and performance of the method: wetting of a chemically heterogeneous substrate in three dimensions, wettingdriven nucleation in a complex two-dimensional domain and three-dimensional diffusion through a porous medium. this is precisely the Young-Dupré angle. When one of the two fluids moves against the other, the contact angle becomes a dynamic quantity, and when the problem is formulated in the framework of conventional hydrodynamics, the contact line motion relative to the solid boundary results in the notorious stress singularity there, as first noted in the pioneering studies by Moffat [43] and Huh and Scriven [34]. Since then there have been numerous analyses and discussions of the singularity over the years, see e.g. [61,32,15] and also recent studies in [58,57] (with the latter one revisiting the classical Cox-Hocking matched asymptotic analysis and providing a correction to it). Recently, it was shown [45] that mesoscopic approaches such as dynamic density functional theory (DDFT) can fix this singularity behavior.A popular model for interface dynamics is the Cahn-Hilliard (CH) equation [10,12], which belongs to the class of phase-field and diffuse interface models. Originally proposed to model spinodal decomposition, the mechanism by which a binary mixture can separate to form two coexisting phases due to, e.g., a change of temperature [12], it has been used in a wide spectrum of different contexts since, such as solidification phenomena [14] and Saffman-Taylor instabilities in Hele-Shaw flows [31]. To account for wetting phenomena and contact lines on solid boundaries, the CH equation can be coupled to a wall boundary condition [11]. Such CH model has been employed successfully in various situations, including microfluidic devices [18, 19, 49, 65], flow in porous media [4], rheological systems [8], and patterning of thin polymer films [38]. Other potential applications include micro-separators [50], fuel cells [3] and CPU chip cooling based on electro-wetting [44]. Many of these applications are characterized by the presence of chemically heterogeneous substrates and/or complex geometries, which make their numerical simulation challenging.The form of the wetting boundary condition is dictated by the form of the wall free energy. For liquid-gas problems linear forms have been adopted, e.g. in the pioneering study by Seppecher [51] and in [9,66]. But a cubic is the lowest-order polynomial required such that the wall free energy can be minimized for the bulk densities, and it prevents the formation of boundary layers on the wall. Cubic forms have been adopted for binary fluid problems, e.g. in [35,69], but also for liquid-gas ones, see [55,56]...