In this paper, we generalize a copula construction method discussed in one of our papers. For this purpose we consider the general form of a linear elliptic PDE. Indeed, a physical interpretation of elliptic equations comes from the notion of conservative flow given by a gradient. This notion provides a mathematical model for equilibrium conservation laws in linear behaviour. This can be applied to many areas of science. Thus, the aim of this paper is to construct a new class of bivariate copulas by solving an elliptic partial differential equation with a Dirichlet condition at the boundary. Copulas belonging to this class allow us to study the stochastic behaviour of the notion of conservative flows. In other words, these copulas will allow us to have an idea on the dependence of those physical phenomena which are governed by elliptic PDEs. For this purpose, we use a discretization method which is the finite difference method which is a common technique for finding approximate solutions of partial differential equations that consists in solving a system of relations (numerical scheme) connecting the values of the unknown functions at some points sufficiently close to each other. For the finite difference method, a mesh is made. This is a set of isolated points called nodes located in the domain of definition of the functions subject to the partial differential equations, a grid on which only the nodes of which the unknowns corresponding to the approximate values of these functions are defined. The mesh also includes nodes located on the boundary of the domain (or at least "close" to this boundary) in order to be able to impose the boundary conditions and/or the initial condition with sufficient accuracy. The primary quality of a mesh is to cover the domain in which it develops as well as possible, to limit the distance between each node and its nearest neighbour. However, the mesh must also allow the discrete formulation of the differentiation operators to be expressed: for this reason, the nodes of the mesh are most often located on a grid whose main directions are the axes of the variables. In the main results of this paper (see section 3), we give a discretization of the solution of the problem followed by a simulation with the MATLAB software of this approximated solution and presenting the discretization errors.