Numerical integration in 2D

I would like to solve the following problem $$\vec{v}(x,y)= k\, \nabla \theta(x,y)$$ with respect to the unknown function $$\theta$$. Parameter $$k$$ is just a real constant quantity.

I have two matrices where I have stored the values of $$v_x(x,y)$$ and $$v_y(x,y)$$, which are the components of the 2D vector field $$\vec{v}$$.

You can convert this equation into a Poisson equation and use standard numerical methods to solve it like finite difference to obtain $$\theta(x,y)$$. If you take a divergence from both side:

$$\nabla^{2} \theta = \frac{1}{k}\nabla \cdot \vec{v}$$

Since, you know the $$\vec{v}$$ you can approximate its divergence as:

$$\nabla \cdot \vec{v} = \frac{\partial v_{x}}{\partial x} + \frac{\partial v_{y}}{\partial y} = \frac{v_{x}(x+\Delta x,y) - v_{x}(x-\Delta x,y)}{2 \Delta x} + \frac{v_{y}(x,y+\Delta y) - v_{y}(x,y-\Delta y)}{2 \Delta y}$$

The only concern here is that how you want to handle the boundaries. If it makes sense for simplest case you can assume a periodic boundary condition here. Now you have $$\nabla \cdot \vec{v}$$ stored in an array and for the rest I assume it is known already. In order to find $$\theta$$ from the defined Poisson equation, you need to know the boundary conditions of $$\theta$$ but again I assume it's defined as periodic boundary conditions on all boundaries. So, now you can discretize the Poisson equation as:

$$\frac{\theta(x+\Delta x,y)+\theta(x-\Delta x,y)-2\theta(x,y)}{\Delta x^{2}} + \frac{\theta(x,y+\Delta y)+\theta(x,y-\Delta y)-2\theta(x,y)}{\Delta y^{2}} = \frac{1}{k} \mathrm{div}(x,y)$$

Where $$\mathrm{div}(x,y)$$ is the array that you stores the values of $$\nabla \cdot \vec{v}$$ in it beforehand. Now you have a standard discretized Poisson equation that you can convert it to a system of linear equations. In fact, I want to convert above system of discretized equations into a matrix form as:

$$A x = b$$

Matrix $$A$$ is called mass matrix and if you have a 2D grid of $$N \times M$$ points, $$A$$'s dimension is $$NM \times NM$$. Now, let's say I can map $$x$$ and $$y$$ grid into matrix coordinate of $$i$$ and $$j$$. In fact, you have $$NM$$ points in your $$x$$ and $$y$$ grid and you can assign an ID from 0 to $$NM - 1$$ to each of them. This would be a function as:

$$F(x,y) = i$$

In fact you want to flatten a 2D grid into 1D as:

$$F(x,y) = x + N y$$

You give $$x$$ and $$y$$ coordinate and it just gives you the id of the point. Then I can assemble $$A$$ and $$b$$ as:

$$A_{F(x,y),F(x,y)} = -2 \Bigg ( \frac{1}{\Delta x^{2}} + \frac{1}{\Delta y^{2}} \Bigg)$$

$$A_{F(x+\Delta x,y),F(x,y)} = \frac{1}{\Delta x^{2}}$$

$$A_{F(x-\Delta x,y),F(x,y)} = \frac{1}{\Delta x^{2}}$$

$$A_{F(x,y + \Delta y),F(x,y)} = \frac{1}{\Delta y^{2}}$$

$$A_{F(x,y - \Delta y),F(x,y)} = \frac{1}{\Delta y^{2}}$$

$$b_{F(x,y)} = \frac{1}{k} \mathrm{div}(x,y)$$

Keep in mind you need to treat boundaries specially based on your boundary condition to assemble $$A$$ and $$b$$. Finally you solve this matrix equation with standard matrix solvers available in Python or any other programming language that you use for your problem.