The Vertex Reinforced Jump Process (VRJP) is a self-reinforcing stochastic process on a graph. There are open questions about its behaviour on Z^d, for example, whether there exists a unique phase transition for d > 2. In order to help build an intuition for these questions, we develop an efficient algorithm for the simulation of the VRJP's random environment on large grid graphs, which can serve as approximations to its dynamics on Z^d. By exploiting the properties of the VRJP's beta-field, we propose an iterative algorithm for its efficient simulation. Its complexity is linear in the number of vertices and cubic in the maximal vertex degree. The random environment can be derived from this beta field by solving a linear system, whose size equals the number of vertices in the graph. To solve this linear system, we use Conjugate Gradients and give a brief discussion of some preconditioning and optimization strategies based on the adjacency matrices of the grid graphs