We propose an idea to solve the Gross–Pitaevskii equation for dark structures inside an infinite constant background density $ρ_∞$=${|ψ_∞|}^2$, without the introduction of artificial boundary conditions. We map the unbounded physical domain $\mathbb{R}^3$ into the bounded domain ${(−1,1)}^3$ and discretize the rescaled equation by equispaced 4th-order finite differences. This results in a free boundary approach, which can be solved in time by the Strang splitting method. The linear part is solved by a new, fast approximation of the action of the matrix exponential at machine precision accuracy, while the nonlinear part can be solved exactly. Numerical results confirm existing ones based on the Fourier pseudospectral method and point out some weaknesses of the latter such as the need of a quite large computational domain, and thus a consequent critical computational effort, in order to provide reliable time evolution of the vortical structures, of their reconnections, and of integral quantities like mass, energy, and momentum. The free boundary approach reproduces them correctly, also in finite subdomains, at low computational cost. We show the versatility of this method by carrying out one- and three-dimensional simulations and by using it also in the case of Bose–Einstein condensates, for which $ψ$→0 as the spatial variables tend to infinity.