This paper considers a family of spatially discrete approximations, including boundary treatment, to initial boundary value problems in evolving bounded domains. The presented method is based on the Cartesian grid embedded Finite-Difference method, which was initially introduced by Abarbanel and Ditkowski [24] [25], for initial boundary value problems on constant irregular domains.We perform a comprehensive theoretical analysis of the numerical issues, which arise when dealing with domains, whose boundaries evolve smoothly in the spatial domain as a function of time. In this class of problems the moving boundaries are impenetrable with either Dirichlet or Neumann boundary conditions, and should not be confused with the class of moving interface problems such as multiple phase flow, solidification, and the stefan problem.Unlike other similar works on this class of problems, the resulting method is not restricted to domains of up to 3-D, can achieve higher then 2 nd -order accuracy both in time and space, and is strictly stable in semi-discrete settings. The strict stability property of the method also implies that the numerical solution remains consistent and valid for a long integration time.A complete convergence analysis is carried in semi-discrete settings, including a detailed analysis for the implementation of the diffusion equation. Numerical solutions of the diffusion equation, using the method for a 2 nd and a 4 th -order of accuracy are carried out in one dimension and two dimensions respectively, which demonstrates the efficacy of the method.