Devising Convergent Numerical Scheme for PDE

I’m currently looking at the PDE
begin{align*}
u_t + left[x(1-y) – (1-x)right]u_x – (1-y) u_y + (z-xy) u_z = (z-xy) u_{xy} – (1-x)u&
end{align*}
with
begin{align*}
u(x,y,z,0) = 1&
(x,y,z) in [0,1]times [0,1] times [0,1]&,
end{align*}
which I found via the conversion of a chemical master equation. The solution $u$ of this PDE is the probability generating function of the master equation. Is there any way to solve this numerically? It resembles advection-diffusion with a source, but it’s not parabolic due to the mixed derivative term and variable coefficients. I don’t believe the lack of boundary conditions is an issue, since if we just look at the advective terms, characteristics flow out of the boundary. Any help/suggestions would be greatly appreciated!