Things get pretty confusing if you attempt to do things in full vectorial notation, so I'll stick exclusively to component notation, with Einstein summations understood
You want to start with the expression
$$
\newcommand{\ue}{\hat{\mathbf e}}\newcommand{\z}{\mathbf 0}
\mathbf A(\mathbf r)=\ue_jA_j(\z)+\ue_jx_k\frac{\partial A_j}{\partial x_k}(\z)+\ue_j x_kx_l\frac{\partial^2 A_j}{\partial x_k\partial x_l}(\z)+ \cdots,
$$
and introduce it into your integral
$$
U=\int \mathbf j(\mathbf r)\cdot\mathbf A(\mathbf r)\mathrm d\mathbf r
$$
to get a multipole series. I will do the monopole and dipole terms and, since you're playing the fuzzy-definition game with the quadrupoles, leave those to you.
The monopole term is easy, and it comes from the zeroth-order term in the integral, which reads
$$
U_0=\int \mathbf j(\mathbf r)\cdot\mathbf A(\z)\mathrm d\mathbf r=\mathbf A(\z)\cdot\int \mathbf j(\mathbf r)\mathrm d\mathbf r,
$$
and which vanishes because you're in a static situation. To prove that this vanishes, consider the $k$th component of the integral of the current, exactly as in this question, giving you
$$
\int j_k(\mathbf r)\mathrm d\mathbf r
=
\int \ue_k\cdot \mathbf j(\mathbf r)\mathrm d\mathbf r
=
\int \left[\nabla\cdot\left(x_k\mathbf j(\mathbf r)\right)-x_k\nabla\cdot\mathbf j(\mathbf r)\right]\mathrm d\mathbf r;
$$
here the second term vanishes because $\nabla\cdot\mathbf j(\mathbf r)=0$, and the first term goes over into the surface integral $\oint_\infty x_k\mathbf j(\mathbf r) \cdot\mathrm d\mathbf S$ at infinity, through the divergence theorem, and this surface integral vanishes since your current is spatially bounded.
The first nontrivial term is the dipole term, which comes from
$$
U_1
= \int j_k(\mathbf r) x_l\frac{\partial A_k}{\partial x_l}(\z) \mathrm d \mathbf r
= \frac{\partial A_k}{\partial x_l}(\z) \int x_lj_k(\mathbf r) \mathrm d \mathbf r,
$$
and this already separates into a field times a moment of the current, but it's obviously not in the form we want it to be in; instead, we want it in the form $\mathbf B(\z)\cdot \mathbf m$, where $\mathbf B = \nabla\times\mathbf A$ and
$$
\mathbf m = \frac12 \int \mathbf r\times\mathbf j(\mathbf r)\mathrm d\mathbf r.
$$
If you work out what we want to get, in all its glory, it's actually not that different, because it reads
\begin{align}
\mathbf B\cdot\mathbf m
& =
\varepsilon_{ikl}\frac{\partial A_l}{\partial x_k}(\z) \,\frac12 \int \varepsilon_{imn}x_mj_n(\mathbf r)\mathrm d\mathbf r
\\& =
\frac12(\delta_{km}\delta_{ln}-\delta_{kn}\delta_{lm})\frac{\partial A_l}{\partial x_k}(\z) \, \int x_mj_n(\mathbf r)\mathrm d\mathbf r
\\& =
\frac12\left[\frac{\partial A_l}{\partial x_k}(\z) \, \int x_kj_l(\mathbf r)\mathrm d\mathbf r
-\frac{\partial A_l}{\partial x_k}(\z) \, \int x_lj_k(\mathbf r)\mathrm d\mathbf r\right],
\end{align}
and this is pretty close to what we have, except for that pesky term with the switched indices. To get a term of that form, you can use an integration-by-parts term as above, which looks as follows: we pull out of our hat the quantity $x_lx_k \,\mathbf j(\mathbf r)$ (much as above we used $x_k \,\mathbf j(\mathbf r)$) and we calculate the integral of its divergence, giving
\begin{align}
\int \nabla\cdot\left(x_lx_k \,\mathbf j(\mathbf r) \right)\mathrm d \mathbf r
& =
\int x_lj_k(\mathbf r) \mathrm d \mathbf r
+\int x_kj_l(\mathbf r) \mathrm d \mathbf r
+\int x_lx_k\,\nabla\cdot\mathbf j(\mathbf r) \mathrm d \mathbf r
.
\end{align}
Here the left-hand side vanishes, because we can pull it to a surface integral at the surface, away from the domain of the current, and the third term on the right also vanishes because the current is static and conserves charge; this means, then, that
$$
\int x_lj_k(\mathbf r) \mathrm d \mathbf r
+\int x_kj_l(\mathbf r) \mathrm d \mathbf r
=0,
$$
and we're essentially done - all that's left to do is to connect the dots.
And, as you can see, trying to do this one order above has nothing but a promise of pain in it, particularly because magnetic quadrupole fields and moments are used relatively rarely and that means that there aren't particularly solid conventions about how exactly one should define the moments, which then makes people a bit more reluctant to use the fields, and it all snowballs down.
To be frank, the magnetic-dipole calculation I've just done (and the electrostatic quadrupole you did in class) is really the highest order to which it makes sense to do the calculation explicitly; if you want to go higher, you really should be setting up a full calculation with arbitrary $l$, with consistent and clear definitions of both the current moments and the external fields at arbitrary multipolarities. Even there, though, it's remarkable that even Jackson forgoes going into those waters, which should tell you something about how worthwhile it is to go down that road.
