You have two options here, neither of which is particularly nice, but whether one is more useful than the other depends on what you have available.
The usual method is to use Fourier series. Working on the square $[0,1] \times [0,1]$ makes the expressions nicer. The boundary conditions say that $v$ is zero on the boundary, so it is appropriate to look for an expansion
$$ v(x,y) = \sum_{n,m=1}^{\infty} a_{n,m} \sin{n\pi x}\sin{m\pi y}. $$
Differentiating, it is easy to see that we then have
$$ -\nabla^2 v = \sum_{n,m=1}^{\infty} (n^2+m^2)\pi^2 a_{n,m} \sin{n\pi x}\sin{m\pi y}. $$
On the other hand,
$$ \int_0^1 \sin{k\pi x}\sin{\ell \pi x} \, dx = \begin{cases} 1/2 & k=\ell \\ 0 & k \neq \ell \end{cases}, $$
so if we multiply the equation by $4\sin{N\pi x}\sin{N\pi y}$ and integrate over the square, we get an expression for $a_{N,M}$:
$$ \int_0^1\int_0^1 (4\sin{N\pi x}\sin{N\pi y})(-\nabla^2v) \, dx \, dy = \dotsb = (N^2+M^2)\pi^2 a_{N,M}. $$
But this must be equal to
$$ \int_0^1\int_0^1 (4\sin{N\pi x}\sin{N\pi y})(c) \, dx \, dy = \frac{(1-(-1)^N)(1-(-1)^M)}{NM\pi^2}c. $$
Then dividing, we find
$$ a_{N,M} = \frac{(1-(-1)^N)(1-(-1)^M)}{NM(N^2+M^2)\pi^4}c $$
and so
$$ v(x,y) = \frac{c}{\pi^4}\sum_{n,m=1}^{\infty} \frac{(1-(-1)^n)(1-(-1)^m)}{nm(n^2+m^2)} \sin{n\pi x}\sin{m\pi y} $$
You can then change variables to get from $[0,1] \times [0,1]$ to $[-B,B] \times [-B,B]$. You'll find that $v$ needs to be rescaled so that the derivative is what it should be.
This is perfectly practical, since you only need a few terms for a very accurate estimate: drawing the graph makes this clear.
An alternative approach would be to use the Dirichlet Green's function for a square of side $a$, which is given by
$$ G(z;w) = \frac{1}{2\pi} \log{\left|\frac{\operatorname{sn}^2{(z/a,1/\sqrt{2})}-\operatorname{sn}^2{(\bar{w}/a,1/\sqrt{2})}}{\operatorname{sn}^2{(z/a,1/\sqrt{2})}-\operatorname{sn}^2{(w/a,1/\sqrt{2})}}\right|}, $$
where sn is a Jacobi elliptic function (see e.g. here for how this comes about). One then has to integrate this over the square to find the solution to Poisson's equation with a constant; as far as I am aware, there's no closed form for this, so you'd have to do it numerically. Theta functions make computing sn easy, but I'm not convinced it's going to be more efficient for your purposes.
