You could try using the Helmholtz decomposition.
If $F$ is a twice-differentiable vector field on a bounded volume $V$ with boundary $S$, then it can be decomposed into divergence-free and curl free components.
$$F=-\nabla\Phi+\nabla\times \mathbf{A}$$
where
$$\Phi(\mathbf{r})=\frac{1}{4\pi}\int_V\frac{\nabla'\cdot F(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}dV'-\frac{1}{4\pi}\oint_S\mathbf{\hat{n}}'\cdot\frac{F(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}dS'$$
$$\mathbf{A}(\mathbf{r})=\frac{1}{4\pi}\int_V\frac{\nabla'\times F(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}dV'-\frac{1}{4\pi}\oint_S\mathbf{\hat{n}}'\times\frac{F(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}dS'$$
$\mathbf{\hat{n}}'$ is the unit outward normal and $\nabla'$ is the gradient with respect to $\mathbf{r}'$ rather than $\mathbf{r}$.
The curl on its own does not have a uniquely-defined inverse. However, the curl and divergence can be combined into a single operator that does have a unique inverse up to boundary conditions.
If the fields are assumed to approach zero at infinity, the boundary integral in each of the above expressions becomes zero. If the field is required to be solenoidal, then the divergence in the first expression will be zero.