I eventually found an answer myself. A function forcing(x,y) has to be written and whatever variable containing the numerical values has to be set as a global variable. Then the MxN numerical matrix has to be transformed into a NumberofNodes-1 array where each line is the value (or interpolated value) of the numerical data at each node.
HOWEVER : I still have difficulties to pass an array of numerical values as Neumann Boundary conditions.. any ideas ?