In the talk we discuss a recently introduced finite element numerical method for the solution of partial differential equations on evolving domains. The approach uses a completely Eulerian description of the domain motion. The physical domain is embedded in a triangulated computational domain and can overlap the time-independent background mesh in an arbitrary way. The numerical method is based on finite difference discretizations of time derivatives and a standard geometrically unfitted finite element method with an additional stabilization term in the spatial domain. The performance and analysis of the method rely on the fundamental extension result in Sobolev spaces for functions defined on bounded domains. Theoretical findings include a complete stability and error analysis, which accounts for discretization errors resulting from finite difference and finite element approximations as well as for geometric errors coming from a possible approximate recovery of the physical domain. We show numerical examples that illustrate the theory and demonstrate the practical efficiency of the method.