All iterative solvers seem to fail for 3D poisson eqn

I have written a function that sets up a sparse matrix A and RHS b for the 3D Poisson equation in a relatively efficient way. The set-up is nothing fancy: I have extended the 2D 5-point stencil to an equivalent 7-point stencil for 3D. Conductors are (at this moment) simply blocks of Dirichlet BCs and I am not (yet) taking dielectrics into account. Therefore, the constructed A is quite straightforward:

• Maximum number of non-zero terms is Nx*Ny*Nz*7
• In practice, my actual number of non-zero terms is about 90% of the maximum

The correctness of my arrays has been verified by running some very simple problems with Nx, Ny and Nz between 50-100 and then using mldivide. However, my peak memory usage is already hitting ~20GB for these test problems, while I'll be using significantly larger arrays for my actual problem.

So now I'm looking into iterative approaches, but my knowledge currently does not extend beyond Gauss-Seidel, Jacobi and SOR. When I attempt to use any of the built-in iterative solvers out of the box, I typically end up with the following message:

METHOD stopped at iteration 2 without converging to the desired tolerance 1e-06
because a scalar quantity became too small or too large to continue computing.
The iterate returned (number 1) has relative residual 0.97.

As this is quite a common system to solve, I was expecting to easily find some resources on this, but I am finding significantly less useful information than I expected. Some advice and pointers would be appreciated!