The isometric embedding problem for Riemannian manifolds, which connects intrinsic and extrinsic geometry, is a central question in differential geometry with deep theoretical significance and wide-ranging applications. Despite extensive analytical progress, the nonlinear and degenerate nature of this problem has hindered the development of rigorous numerical analysis in this area. As the first step toward addressing this gap, we study the numerical approximation of Weyl’s problem, i.e., the isometric embedding of two-dimensional Riemannian manifolds with positive Gaussian curvature into $\mathbb{R}^3$, by establishing a new weak formulation that naturally leads to a numerical scheme well suited for high-order finite element discretization, and conducting a systematic analysis to prove the well-posedness of this weak formulation, the existence and uniqueness of its numerical solution, as well as its convergence with error estimates. This provides a foundational framework for computing isometric embeddings of Riemannian manifolds into Euclidean space. Numerical experiments support the theoretical analysis by demonstrating the convergence of the method and its effectiveness in simulating isometric embeddings of given Riemannian manifolds.