I think you want to solve Poisson's equation in a rectangle with Neumann boundary conditions  on two sides or all four sides. Any way apply 5-point finite difference scheme and form a block tri-diagonal system and solve it by Gauss Seidel method. Since it is not time dependent problem, so nothing like explicit method. For detail see, G D Smith or Greenspan or L Collatz or M K Jain book.

Use parallel computing to reduce the time. It works well for this type of problems.

see also FEM in the following link with c++ codes:

https://people.sc.fsu.edu/~jburkardt/cpp_src/fem2d_poisson_sparse/fem2d_poisson_sparse.html

Lots of thanks to all colleagues!

Dear R.C.Mittal your advice is right what I was looking for.

Dear Seyed Hamed Hashemi Mehne, thanks a lot, I was planing to go parallel my whole code when it will be working ok!

Since after discretisation you will got a linear system to solve A*\phi = f,  where A is symetric positive semidefinite and f is the rhsm I think you still can use preconditioned conjugate gradient method  (for instance CG-ILU) .

Best

It should be noted that solutions to this problem are not unique. Take any solution and add a constant and you should get a new valid solution.  That being said some numerical schemes will have difficulty with this.  I would suggest a FEM with some additional constraint that the integral of the solution over the domain is some constant like 0.

Poisson equation is an elliptic equation and hence it strongly depends on the boundary condition. As your domain of solution is of rectangular shape, therefore it is easy to use Neumann  boundary condition on two parallel sides and Dirichlet boundary condition on other two sides. Since the solution of poisson equation is obtained through iterative process, therefore it is easy to implement successive over relaxation (SOR). SOR quiet easy to use and faster than GS and Jacobi method.