Natural neighbor coordinates [20] are optimum weighted-average measures for an irregular arrangement of nodes in R n . [26] used the notion of BΓ©zier simplices in natural neighbor coordinates Ξ¦ to propose a C 1 interpolant. The C 1 interpolant has quadratic precision in β¦ β R 2 , and reduces to a cubic polynomial between adjacent nodes on the boundary ββ¦. We present the C 1 formulation and propose a computational methodology for its numerical implementation (Natural Element Method) for the solution of partial differential equations (PDEs). The approach involves the transformation of the original Bernstein basis functions B 3 i (Ξ¦) to new shape functions Ξ¨(Ξ¦), such that the shape functions Ο3I-2(Ξ¦), Ο3I-1(Ξ¦), and Ο3I (Ξ¦) for node I are directly associated with the three nodal degrees of freedom wI , ΞΈI x , and ΞΈI y , respectively. The C 1 shape functions interpolate to nodal function and nodal gradient values, which renders the interpolant amenable to application in a Galerkin scheme for the solution of fourth-order elliptic PDEs. Results for the biharmonic equation with Dirichlet boundary conditions are presented. The generalized eigenproblem is studied to establish the ellipticity of the discrete biharmonic operator, and consequently the stability of the numerical method.