A general multilayer diffusion problem subject to inhomogeneous boundary conditions with time-dependent boundaries is considered. The embedding method is used to find the analytical solution, which is expressed in terms of time-dependent functions that satisfy a system of linear Volterra integral equations of the first kind. A boundary element method is developed to numerically solve the integral equations and results of numerical simulations are presented. The proposed approach is applicable to multilayer diffusion problems defined on bounded and unbounded spatial domains.