In this paper, we investigate inhomogeneous and simultaneous Diophantine approximation in beta dynamical systems. For $\beta>1$ let $T_{\beta}$ be the $\beta$-transformation on $[0,1]$. We determine the Lebesgue measure and Hausdorff dimension of the set \[\left\{(x,y)\in [0,1]^2: |T_{\beta}^nx-f(x,y)|<\varphi(n)\text{ for infinitely many }n\in\mathbb{N}\right\},\] where $f:[0,1]^2\to [0,1]$ is a Lipschitz function and $\varphi$ is a positive function on $\mathbb{N}$. Let $\beta_2\geq \beta_1>1$, $f_1,f_2:[0,1]\to [0,1]$ be two Lipschitz functions, $\tau_1,\tau_2$ be two positive continuous functions on $[0,1]$. We also determine the Hausdorff dimension of the set \[\left\{(x,y)\in [0,1]^2: \begin{aligned}&|T_{\beta_1}^nx-f_1(x)|<\beta_1^{-n\tau_1(x)}\\ &|T_{\beta_2}^ny-f_2(y)|<\beta_2^{-n\tau_2(y)}\end{aligned}\text{ for infinitely many }n\in\mathbb{N}\right\}.\] Under certain additional assumptions, the Hausdorff dimension of the set \[\left\{(x,y)\in [0,1]^2: \begin{aligned}&|T_{\beta_1}^nx-g_1(x,y)|<\beta_1^{-n\tau_1(x)}\\ &|T_{\beta_2}^ny-g_2(x,y)|<\beta_2^{-n\tau_2(y)}\end{aligned}\text{ for infinitely many }n\in\mathbb{N}\right\}\] is also determined, where $g_1,g_2:[0,1]^2\to [0,1]$ are two Lipschitz functions.