We consider a $d$-dimensional random field $u=(u(x), x\in D)$ that solves a system of elliptic stochastic equations on a bounded domain $D\subset \mathbb{R}^k$, with additive white noise and spatial dimension $k=1,2,3$. Properties of $u$ and its probability law are proved. For Gaussian solutions, using results from [Dalang and Sanz-Solé, 2009], we establish upper and lower bounds on hitting probabilities in terms of the Hausdorff measure and Bessel-Riesz capacity, respectively. This relies on precise estimates on the canonical distance of the process or, equivalently, on $L^2$ estimates of increments of the Green function of the Laplace equation.