The diffuse-domain, or smoothed boundary, method is an attractive approach for solving partial differential equations in complex geometries because of its simplicity and flexibility. In this method the complex geometry is embedded into a larger, regular domain. The original PDE is reformulated using a smoothed characteristic function of the complex domain, here constructed from a level-set (signed distance) function, and source terms are introduced to approximate the boundary conditions. The reformulated equation can be solved by standard numerical methods and the same solver can be used for any domain geometry. A challenge is making the method higher order accurate. For Dirichlet boundary conditions, which we focus on here, current implementations demonstrate a wide range in their accuracy but generally the methods yield at best first order accuracy in $\epsilon$, the parameter that characterizes the width of the region over which the characteristic function is smoothed. Typically, $\epsilon\propto h$, the grid size. Here, we analyze the diffuse-domain PDEs using matched asymptotic expansions and explain the observed behaviors. Our analysis also identifies simple modifications to the diffuse-domain PDEs that yield higher-order accuracy in $\epsilon$, e.g., $O(\epsilon^2)$ in the $$L^2$$ norm and $O(\epsilon^p)$ with $1.5\le p\le 2$ in the $$L^{\infty}$$ norm. Our analytic results are confirmed numerically in stationary and moving domains. Finally, we extend our results to two-phase fluid problems and propose a diffuse domain approach to simulate two-phase flows in complex geometries.