Hi! There are multiple packages than can do multivariate interpolations. You can have a look at the following thread: 3D interpolation in Julia. To name some you might want to look at: Dierckx.jl, GridInterpolations.jl, GMT.jl, DIVAnd.jl, RadialBasisFunctions.jl, ScatteredInterpolation.jl, KernelInterpolation.jl. Which of these is suited the best for your use case depends on how your data looks like, especially the location of your data points in \mathbb{R}^2 (e.g., grid data or scattered data). If you have any questions about one (or multiple) of these packages, feel free to ask.

# x coord; dim = n x 1
X = [0.0 1.0 3.0 0.0 1.0 3.0 0.0 1.0 3.0]
# y coord; dim = n x 1
Y = [0.0 0.0 0.0 1.0 1.0 1.0 3.0 3.0 3.0]
# PDE sol; dim = n x 1
sol = [0.5 1.0 1.5 1.0 0.5 1.0 1.5. 2.0 0.5]

(X[i], Y[i]) is a node in a grid and sol[i] is the solution at (X[i], Y[i]).
I want to build now an interpolation of the vector sol on the whole domain.

AFAIU for the answers and the other threads, I need a a n x n matrix for sol?

If your data is from a grid, you will get the most efficiency from providing sol as a matrix that conforms to the grid. Do you have data on the entire grid, or do you need a scattered interpolant?

If you do actually want to interpolate with data like from your example, you can do e.g.

julia> using NaturalNeighbours
julia> X = [0.0 1.0 3.0 0.0 1.0 3.0 0.0 1.0 3.0] |> vec;
julia> Y = [0.0 0.0 0.0 1.0 1.0 1.0 3.0 3.0 3.0] |> vec;
julia> sol = [0.5 1.0 1.5 1.0 0.5 1.0 1.5 2.0 0.5] |> vec;
julia> itp = interpolate(X, Y, sol); # also see ?interpolate for some keywords
julia> itp(0.2, 0.3)
0.69
julia> itp([0.5, 0.2, 0.3], [0.1, 0.2, 0.9]; method = Laplace()) # see ?interpolate
3-element Vector{Float64}:
0.75
0.66
0.8300000000000001

If you are trying to interpolate a smooth function f(x,y), e.g. the solution of a PDE as a function of some parameters, and you have the freedom to choose the x,y points, then you will probably be much better off choosing the x,y points to lie on a Chebyshev grid and then use Chebyshev interpolation, which converges exponentially fast for smooth functions. See FastChebInterp.jl.

On the other hand if your x,y points come from how you solve the PDE, e.g. a finite-difference grid, then you should choose your interpolation based on what computational method you used. e.g. if you used 2nd-order centered differences then you should just use bilinear interpolation, and if you used finite elements you should use the FEM basis functions.