Statistical mechanics
Van der Waals equation can be derived from virial/cluster expansion in terms of the average/mean density of the liquid/gas:
$$
\rho(\mathbf{r})=\left\langle \sum_{i=1}^N\delta(\mathbf{r}-\mathbf{r}_i)\right\rangle,
$$
where $\mathbf{r}_i$ is the position of the i-th atom/molecule.
Thermodynamics
Another way to look at it is that pressure is the derivative of the free energy in respect to volume, so one can obtain the expression for the free energy by integrating the pressure expressed from van der Waals equation in respect to volume. One can then analyze liquid-gas transition in terms of Landau theory of phase transitions, with density serving as the order parameter.
See, e.g., section 8.5 in Lecture notes on Statistical Physics by Lydéric Bocquet
See also Wikipedia for the expression for the Helmholz free energy in more conventional terms.
Simple mean-field derivation of Van der Waals equation
Simple mean-field derivation of VdW equation can be performed along the lines of the derivation given in Wikipedia. I sketch it here in even more simplified/clarified form.
The energy of a gas of N interacting particles is
$$
H(\mathbf{x}_1, \mathbf{p}_1;...;\mathbf{x}_N, \mathbf{p}_N)=
\sum_{i=1}^N\frac{\mathbf{p}_i^2}{2m}+\frac{1}{2}\sum_{i=1}^N\sum_{j=1,j\neq i}^NU(|\mathbf{x}_i-\mathbf{x}_j|),
$$
where the potential is attractive, i.e., $U(|\mathbf{x}_i-\mathbf{x}_j|)<0$ (unless the molecules/atoms are very close to each other)
The partition function is
$$
Z=\int \prod_{i=1}^Nd^3\mathbf{x}_id^3\mathbf{p}_ie^{-\beta H(\mathbf{x}_1, \mathbf{p}_1;...;\mathbf{x}_N, \mathbf{p}_N)}=
\int d^3\mathbf{p}_1...d^3\mathbf{p}_Ne^{-\beta \sum_{i=1}^N\frac{\mathbf{p}_i^2}{2m}}Z_c=\\\left(2\pi mk_B T\right)^{\frac{3N}{2}}Z_c
$$
where the configuration integral is
$$
Z_c=\int d^3\mathbf{x}_1...d^3\mathbf{x}_Ne^{-\frac{\beta}{2}\sum_{i=1}^N\sum_{j=1,j\neq i}^NU(|\mathbf{x}_i-\mathbf{x}_j|)},
$$
and $\beta=1/(k_BT)$.
We now perform the mean-field approximation replacing the pairwise interactions between particles by an interaction of each particle with an averaged field of all other particles:
$$
\frac{1}{2}\sum_{i=1}^N\sum_{j=1,j\neq i}^NU(|\mathbf{x}_i-\mathbf{x}_j|)\approx
\frac{1}{2}\sum_{i=1}^N \frac{1}{V}\int_V d^3\mathbf{x_j}U(|\mathbf{x}_i-\mathbf{x}_j|)\approx\\
\frac{1}{2}\sum_{i=1}^N \frac{N}{V}\int_V d^3\mathbf{y}U(|\mathbf{y}|)=-\frac{a'N^2}{V},
$$
where the minus sign is due to the fact that the potential is attractive, so that $a'>0$.
The configuration integral is then
$$
Z_c=\int d^3\mathbf{x}_1...d^3\mathbf{x}_N e^{\beta\frac{a'N^2}{V}}=
(V-b'N)^N e^{\beta\frac{a'N^2}{V}},
$$
where we took account for the fact that the particles cannot be too close to each other, and hence the volume integral for each particle is reduced by the factor accounting for the volume of all other particles.
We now can write down the free energy of the gas as
$$
F=-\frac{1}{\beta}\log Z=-\frac{a'N^2}{V}-\frac{N}{\beta}\log(V-b'N)-\frac{3N}{2\beta}\log\left(2\pi mk_BT\right)
$$
The gas pressure is given by (since $dF=-SdT-pdV$)
$$
P=-\left(\frac{\partial F}{\partial V}\right)_T=
\frac{Nk_BT}{V-Nb'}-\frac{a'N^2}{V^2}
$$
Converting to molar fractions $n=N/N_A$, $b=N_Ab'$, $a=a'N_A^2$, we obtain the usual form of the VdW equation
$$
P=\frac{nRT}{V-nb}-\frac{an^2}{V^2}\Leftrightarrow
\left(P+\frac{an^2}{V^2}\right)(V-nb)=nRT.
$$
Best Answer
This feels a bit like a homework question, so you will only get half of the answer and need to do the final part. Make sure you understand what is going on though, so you can finish fast!
You need to try and write the VdW equation in the form suggested, and then find terms $\sim n/V$ or $\sim n^2 / V^2$ and identify the coefficients ($B_1$ and $B_2$ in this case).
You can start by rewriting
$$(p+{an^2 \over V^2})(V-nb) = nRT$$ as $$pV-pnb+{an^2\over V}-{abn^3\over V^2} = nRT$$
so that
$$p(V-nb)=nRT+{an^2\over V}+{abn^3\over V^2}$$
Notice that we had to keep the $p$-dependet terms on the left, because pressure has to only appear on the left. On the other hand, $V$ can appear on the right but only in powers of $n/V$. So we need to decouple things.
To do that, we first collect all things on the right:
$$p(V-nb)=nRT\left(1+{a \over RT} {n\over V}+{ab\over RT}{n^2\over V^2}\right)$$
and you see we only get powers of $n/V$ inside the parenthesis.
We would be done if we did not have that annoying $p(V-nb)$ term that is keeping an extra $..-nb$ on the left.
To get rid of that, we can use the fact that
$$p(V-nb)=pV(1-b{n\over V})$$
and we can now isolate $pV$ and divide by the remaining term
$$pV = nRT {1\over 1-b{n\over V}}\left(1+{a \over RT} {n\over V}+{ab\over RT}{n^2\over V^2}\right)$$
And you see that we have $pV$ on the left (good) and only powers of $n/V$ on the right (good). Unfortunately, one $n/V$ is at the denominator still...
So all that is left - and I leave it to you so you can finish the exercise - is to expand ${1\over 1-b{n\over V}}$ in series (because we assume $nb\ll V$ of course, as it is a correction) and again collect terms in order $\sim n / V$, $\sim n^2 / V^2$, $\sim n^3 / V^3$ etc.