We present a natural element method to treat higher-order spatial derivatives in the Cahn-Hilliard equation. The Cahn-Hilliard equation is a fourth-order nonlinear partial differential equation that allows to model phase separation in binary mixtures. Standard classical C0-continuous finite element solutions are not suitable because primal variational formulations of fourth-order operators are well-defined and integrable only if the finite element basis functions are piecewise smooth and globally C 1-continuous. To ensure C1-continuity, we develop a natural-element-based spatial discretization scheme. The C1- continuous natural element shape functions are achieved by a transformation of the classical Farin interpolant, which is basically obtained by embedding Sibsons natural element coordinates in a Bernstein-Bézier surface representation of a cubic simplex. For the temporal discretization, we apply the (second-order accurate) trapezoidal time integration scheme supplemented with an adaptively adjustable time step size. Numerical examples are presented to demonstrate the efficiency of the computational algorithm in two dimensions. Both periodic Dirichlet and homogeneous Neumann boundary conditions are applied. Also constant and degeneratemobilities are con-sidered.We demonstrate that the use of C1-continuous natural element shape functions enables the computation of topologically correct solutions on arbitrarily shaped domains.