Including even one Dirichelt condition changes the problem you are trying to solve and will not give you the correct solution! You must fulfil the Discrete Compatibility Criteria, see e.g. first and second pages of:
It basically states that the summation of each element of b must be equal to zero for all-Neumann problems, otherwise your solution will drift as there is nothing holding it at the boundaries.
You simply need to add a line stating that b = b - mean(b) in your code before solving, if you are solving in double precision. In single precision it might also be necessary to ensure, at every other iteration, that the residual in biCGstab also meets the Discrete Compatibility Criteria.
I attempted to solve an all-Neumann problem using CG by imposing a Dirichlet condition at a single point, as one is usually told to do anecdotally by e.g. some Prof., and then compared it to using a Discrete Cosine Transform approach. It does not yield the same result, using the Discrete Compatibility Criteria instead does.