diff --git a/docs/make.jl b/docs/make.jl
index 602c607..74afab5 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -18,6 +18,7 @@ makedocs(
         "Examples" => [
             "examples/normal_distribution.md",
             "examples/general_distribution.md",
+            "examples/multivariate_distribution.md",
         ],
         "API Reference" => "api_reference.md",
     ]
diff --git a/docs/src/examples/multivariate_distribution.md b/docs/src/examples/multivariate_distribution.md
new file mode 100644
index 0000000..03a5c89
--- /dev/null
+++ b/docs/src/examples/multivariate_distribution.md
@@ -0,0 +1,81 @@
+# Example: Multivariate Distribution
+
+This example demonstrates the use of `PointEstimateMethod.jl` to approximate a set of independent univariate distributions jointly using the multivariate PEM.
+
+## Setup
+
+```julia
+using Distributions
+using PointEstimateMethod
+```
+
+## Approximating Two Independent Distributions
+
+We consider two independent distributions:
+- A Normal distribution: ``X_1 \sim \mathcal{N}(2.0, 0.5)``
+- A LogNormal distribution: ``X_2 \sim \text{LogNormal}(0.5, 0.3)``
+
+```julia
+d1 = Normal(2.0, 0.5)
+d2 = LogNormal(0.5, 0.3)
+
+distributions = [d1, d2]
+N_pem = 3  # Number of point estimate points per distribution
+
+pem_output = pem(distributions, N_pem)
+```
+
+The output `pem_output` is a named tuple with two fields:
+
+- `pem_output.x` — a `K × N` matrix of location points (one row per distribution)
+- `pem_output.p` — a `K × N` matrix of probability weights (one row per distribution)
+
+```julia
+println("Locations matrix (K × N):")
+println(pem_output.x)
+
+println("Probabilities matrix (K × N):")
+println(pem_output.p)
+```
+
+## Verifying Moment Preservation
+
+The sum of probabilities for each distribution should equal `1/K`:
+
+```julia
+K = length(distributions)
+for k in 1:K
+    println("Sum of probabilities for distribution $k: ", sum(pem_output.p[k, :]))  # ≈ 1/K
+end
+```
+
+The weighted mean for each distribution should match the true mean:
+
+```julia
+for k in 1:K
+    weighted_mean = sum(pem_output.p[k, j] * pem_output.x[k, j] for j in 1:N_pem)
+    println("Distribution $k — PEM weighted mean: $(real(weighted_mean)), true mean: $(mean(distributions[k]))")
+end
+```
+
+## Using Sample Data
+
+The multivariate PEM also accepts a `K × M` matrix of samples, where `K` is the number of variables and `M` is the number of samples:
+
+```julia
+using Random
+Random.seed!(42)
+
+M = 10000
+samples = vcat(
+    rand(Normal(2.0, 0.5), 1, M),
+    rand(LogNormal(0.5, 0.3), 1, M)
+)  # 2 × 10000 matrix
+
+pem_samples = pem(samples, 3)
+
+println("Locations from samples:")
+println(pem_samples.x)
+println("Probabilities from samples:")
+println(pem_samples.p)
+```
diff --git a/docs/src/index.md b/docs/src/index.md
index c73352c..edb47ce 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -37,6 +37,6 @@ pem_output.p  # probability of each point
 ## Contents
 
 ```@contents
-Pages = ["index.md", "installation.md", "examples/normal_distribution.md", "examples/general_distribution.md", "api_reference.md"]
+Pages = ["index.md", "installation.md", "examples/normal_distribution.md", "examples/general_distribution.md", "examples/multivariate_distribution.md", "api_reference.md"]
 Depth = 2
 ```
diff --git a/src/PointEstimateMethod.jl b/src/PointEstimateMethod.jl
index 7e9ef9d..8e6009a 100644
--- a/src/PointEstimateMethod.jl
+++ b/src/PointEstimateMethod.jl
@@ -23,6 +23,7 @@ DEFAULT_SOLVER = optimizer_with_attributes(
 )
 
 include("auxiliaries.jl")
-include("pem.jl")
+include("pem_univariate.jl")
+include("pem_multivariate.jl")
 
 end
diff --git a/src/pem_multivariate.jl b/src/pem_multivariate.jl
new file mode 100644
index 0000000..4fa0ebc
--- /dev/null
+++ b/src/pem_multivariate.jl
@@ -0,0 +1,242 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2022 Davide Fioriti
+#
+# SPDX-License-Identifier: GPL-3.0-or-later
+# coding: utf-8
+
+
+"""
+    pem(d, N; mean_fun=mean, central_moment_fun=moment, optimizer=HiGHS.Optimizer)
+
+Point Estimate Method to identify N estimate points for the univariate distribution d.
+
+Parameters
+----------
+- d :: UnivariateDistribution
+    Distribution under interest
+- N :: Integer
+    Number of desired estimate points
+- mean_fun :: Function (optional)
+    Function used to calculate the mean value of the distribution
+- central_moment_fun :: Function (optional)
+    Function used to calculate the central moment of the distribution d
+- montecarlo_sampling :: Integer (optional, default 1e6)
+    Number of Monte Carlo samples used in the sampling procedure if a non-specific moment function is available
+- optimizer (optional)
+    JuMP optimizer for executing the optimization
+
+Returns
+-------
+- (x, p) :: NamedTuple
+    - x :: Vector
+        Return the location points
+    - p :: Vector
+        Return the probability of each point
+
+"""
+function pem(
+        d::Vector{<:UnivariateDistribution},
+        N::Integer;
+        mean_fun::Function=Distributions.mean,
+        central_moment_fun::Function=Distributions.moment,
+        montecarlo_sampling::Integer= 1000000,
+        optimizer=DEFAULT_SOLVER,
+    )
+
+    K = length(d)
+    
+    if all(hasmethod(central_moment_fun, Tuple{typeof(d[k]), Int}) for k in 1:K)
+    
+        # central moments
+        m_list = Dict(
+            (k,i)=>central_moment_fun(d[k], i)
+            for k = 1:K for i = 1:(2*N)
+        )
+
+        mean_values = [mean_fun(d[k]) for k in 1:K]
+
+        return pem(mean_values, m_list, N; optimizer=optimizer)
+    else
+        @info """Function $(string(central_moment_fun)) does not have a direct implementation for Distribution $(string(d)). Perform Monte Carlo sempling over the distribution with $montecarlo_sampling points"""
+        sampled_set = zeros(K, montecarlo_sampling)
+        for k = 1:K
+            sampled_set[k, :] = rand(d[k], montecarlo_sampling)
+        end
+
+        return pem(sampled_set, N; optimizer=optimizer)
+    end
+end
+
+
+"""
+    pem(d, N; mean_fun=mean, central_moment_fun=moment, optimizer=HiGHS.Optimizer)
+
+Point Estimate Method to identify N estimate points for an experimental multi-variate distribution
+represented by the array of elements d.
+The first dimension of the array is the number of variables, while the second dimension is the number of samples for each variable.
+
+Parameters
+----------
+- d :: Array{<:Real, 2}
+    Distribution under interest
+- N :: Integer
+    Number of desired estimate points
+- mean_fun :: Function (optional)
+    Function used to calculate the mean value of the distribution
+- central_moment_fun :: Function (optional)
+    Function used to calculate the central moment of the distribution d
+- optimizer (optional)
+    JuMP optimizer for executing the optimization
+
+Returns
+-------
+- (x, p) :: NamedTuple
+    - x :: Array
+        Return the location points
+    - p :: Array
+        Return the probability of each point
+
+"""
+function pem(
+        d::Array{<:Real, 2},
+        N::Integer;
+        mean_fun::Function=Distributions.mean,
+        central_moment_fun::Function=Distributions.moment,
+        optimizer=DEFAULT_SOLVER,
+    )
+    
+    ## Execution
+    ## Solving methodology by https://www.jstor.org/stable/2631060
+
+    K = size(d, 1)
+
+    mean_values = [mean_fun(d[k,:]) for k in 1:K]
+    
+    # moments
+    m_list = Dict(
+        (k, i)=>central_moment_fun(d[k,:], i)
+        for k = 1:K for i = 1:(2*N)
+    )
+
+    return pem(mean_values, m_list, N; optimizer=optimizer)
+end
+
+
+"""
+    pem(mean_values, d, m_list, N; optimizer=HiGHS.Optimizer)
+
+Point Estimate Method to identify estimate points for K independent univariate distributions with moments given by m_list and mean values given by mean_values.
+
+This function is based on the methodology proposed by:
+- H.P.Hong, An efficient point estimate method for probabilistic analysis, Reliability Engineering and System Safety, 1998, https://doi.org/10.1016/S0951-8320(97)00071-9
+- Miller, Allen C., and Thomas R. Rice. “Discrete Approximations of Probability Distributions.” Management Science 29, no. 3 (1983): 352–62. http://www.jstor.org/stable/2631060.
+
+Parameters
+----------
+- mean_value :: Vector{<:Real}
+    Mean value of the distribution
+- m_list :: Dict{Tuple{<:Integer,<:Integer}, <:Real}
+    Dictionary representing the central moments of the distribution.
+    The keys of the dictionary are tuples (k, m) where k is the distribution index and m is the moment order.
+    The value corresponds to the value of the moment.
+    Note: they are central moments referred to the mean. As such, the moment of order 1 is 0.0.
+- N :: Integer
+    Number of desired estimate points
+- optimizer (optional)
+    JuMP optimizer for executing the optimization
+
+Returns
+-------
+- (x, p) :: NamedTuple
+    - x :: Vector
+        Return the location points
+    - p :: Vector
+        Return the probability of each point
+
+"""
+function pem(
+        mean_values::Vector{<:Real},
+        m_list::Dict{<:Tuple{<:Integer,<:Integer}, <:Real},
+        N::Integer;
+        optimizer=DEFAULT_SOLVER,
+    )
+
+    # number of distributions
+    K = length(mean_values)
+    
+    # ensure consistency of the input
+    expected_keys = Set((k, i) for k in 1:K for i in 1:(2*N))
+    @assert Set(keys(m_list)) == expected_keys "The input moment dictionary does not match the expected index in th form (k,m) where k is the distribution index and m is the moment order."
+    
+    ## Execution
+    ## Solving methodology by https://www.jstor.org/stable/2631060
+    
+    # lambda i value
+    # The moment of order 0 has value 1/K, as we are considering K independent distributions
+    λ = Dict(
+        (k,i)=>((i==0) ? 1. / K : m_list[(k,i)])
+        for k=1:K for i = 0:2*N
+    )
+    
+    ## 1) Preliminary model to get the coefficients of polynomial described in section 4
+    ##    of https://www.jstor.org/stable/2631060
+    
+    model = Model(optimizer)
+    
+    # coefficients of auxiliary polynomial \sum_{k=0}^N C_k x^k = π(x) = (x - x_1) ... (x - x_N) for each distribution
+    @variable(model, C[k=1:K,i=0:N-1])
+
+    @constraint(
+        model,
+        aux_poly_balance[k=1:K,i=0:N-1],
+        sum(
+            C[k,p]*λ[k,p+i]
+            for p=0:N-1
+        ) == -λ[k,N+i]
+    )
+    
+    # Determine coefficients
+    optimize!(model)
+    
+    # 2) postprocess the coefficients to obtain the desired locations
+    ϵ = zeros(K, N)
+
+    for k in 1:K
+        # get the coefficients of the polynomial
+        poly_coeffs = [[value(C[k,i]) for i = 0:N-1]; 1.0]
+        # get the roots of the polynomial and get the standardized locations of the distribution
+        poly = Polynomials.Polynomial(poly_coeffs)
+        ϵ[k,:] = Polynomials.roots(poly)
+    end
+    
+    # obtain the probabilities of such locations
+    postmodel = Model(optimizer)
+    
+    @variable(postmodel, probabilities[k=1:K,j=1:N])
+
+    # balance the moments for each distribution
+    @constraint(
+        postmodel,
+        balance_moments[k=1:K,i=1:(2*N-1)],
+        sum(probabilities[k,j] * ϵ[k,j]^i for j in 1:N) == λ[k,i]
+    )
+    
+    # set the probabilities to sum to 1/K for each distribution, as we are considering K independent distributions
+    @constraint(
+        postmodel,
+        balance_probabilities[k=1:K],
+        sum(probabilities[k,:]) == 1. / K
+    )
+
+    # optimize the model
+    optimize!(postmodel)
+    
+    # get probabilities
+    p = reshape(value.(probabilities), K, N)
+
+    # get locations of the estimated points
+    x = ϵ .+ reshape(mean_values, K, 1)
+    
+    # results = NamedTuple{(:x, :p)}.(zip(x, p))
+    return (x=x, p=p)
+end
diff --git a/src/pem.jl b/src/pem_univariate.jl
similarity index 99%
rename from src/pem.jl
rename to src/pem_univariate.jl
index 4d0d33e..6e2628b 100644
--- a/src/pem.jl
+++ b/src/pem_univariate.jl
@@ -44,7 +44,7 @@ function pem(
     
     if hasmethod(central_moment_fun, Tuple{typeof(d), Int})
     
-        # lambda i value
+        # central moments
         m_list = Dict(
             i=>central_moment_fun(d, i)
             for i = 1:(2*N)
@@ -59,6 +59,7 @@ function pem(
     end
 end
 
+
 """
     pem(d, N; mean_fun=mean, central_moment_fun=moment, optimizer=HiGHS.Optimizer)
 
@@ -204,4 +205,4 @@ function pem(
     
     # results = NamedTuple{(:x, :p)}.(zip(x, p))
     return (x=x, p=p)
-end
\ No newline at end of file
+end
diff --git a/test/Examples.jl b/test/Examples.jl
index 07a2839..9dd21a6 100644
--- a/test/Examples.jl
+++ b/test/Examples.jl
@@ -13,6 +13,10 @@ module Examples
         Example("Normal_3", Normal(), 3),
         Example("Normal_truncated_3", truncated(Normal(1.0, 0.4), 0.0, +Inf), 3),
         Example("Normal_truncated_9", truncated(Normal(1.0, 0.4), 0.0, +Inf), 9),
+        Example("MultiVariate_Normal_2", [Normal()], 2),
+        Example("MultiVariate_2Normal_3", [Normal(), Normal()], 3),
+        Example("MultiVariate_3Normal_2", [Normal(), Normal(), Normal()], 2),
+        Example("MultiVariate_3Normal_3", [Normal(), Normal(), Normal()], 3),
     ]
 
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index b56922c..a767e63 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -11,11 +11,18 @@ ATOL_TEST = 1e-4
 RTOL_TEST = 1e-4
 SEED = 0
 
-function vector_approx_test(x_test::Vector, x_validation::Vector, rtol=RTOL_TEST, atol=ATOL_TEST)
+# function to test the approximation of vectors
+function approx_test(x_test::Vector, x_validation::Vector, rtol=RTOL_TEST, atol=ATOL_TEST)
     length(x_test) != length(x_validation) && return false
     return all(isapprox(xt, xv; rtol=rtol, atol=atol) for (xt, xv) in zip(x_test, x_validation))
 end
 
+# function to test the approximation of matrices
+function approx_test(x_test::Matrix, x_validation::Matrix, rtol=RTOL_TEST, atol=ATOL_TEST)
+    size(x_test) != size(x_validation) && return false
+    return all(isapprox(xt, xv; rtol=rtol, atol=atol) for (xt, xv) in zip(x_test, x_validation))
+end
+
 "Function to test examples"
 function test_example(example_name, testing_function, args...)
 
@@ -27,21 +34,42 @@ function test_example(example_name, testing_function, args...)
 
     path_solution = joinpath(BASE_FOLDER, "test", "testcases", string(testing_function), example_name * ".yml")
 
-    I_calc = sortperm(calc_solution.x)
+    
+    if calc_solution.x isa Vector
+        # Univariate distribution
+
+        I_calc = sortperm(calc_solution.x)
+        x_sorted = calc_solution.x[I_calc]
+        p_sorted = calc_solution.p[I_calc]
+
+    else
+        # Multivariate distribution
+        K = size(calc_solution.x, 1)
+        I_calc = [sortperm(calc_solution.x[k, :]) for k in 1:K]
+        x_sorted = Matrix(reduce(hcat, [calc_solution.x[k, I_calc[k]] for k in 1:K])')
+        p_sorted = Matrix(reduce(hcat, [calc_solution.p[k, I_calc[k]] for k in 1:K])')
+    end
     
     if isfile(path_solution)
         # if the file exists run tests
         proven_solution = YAML.load_file(path_solution)
+
+        x_validation = proven_solution["x"] isa Vector{<:Vector} ? Matrix(reduce(hcat, proven_solution["x"])') : proven_solution["x"]
+        p_validation = proven_solution["p"] isa Vector{<:Vector} ? Matrix(reduce(hcat, proven_solution["p"])') : proven_solution["p"]
         
-        @test vector_approx_test(calc_solution.x[I_calc], proven_solution["x"])
-        @test vector_approx_test(calc_solution.p[I_calc], proven_solution["p"])
+        @test approx_test(x_sorted, x_validation)
+        @test approx_test(p_sorted, p_validation)
     else
         # otherwise create the tests
         mkpath(dirname(path_solution))
 
+        x_for_yaml = x_sorted isa Matrix ? [x_sorted[i, :] for i in axes(x_sorted, 1)] : x_sorted
+        p_for_yaml = p_sorted isa Matrix ? [p_sorted[i, :] for i in axes(p_sorted, 1)] : p_sorted
+
+
         dict_calc_solution = Dict(
-            "x"=>calc_solution.x[I_calc],
-            "p"=>calc_solution.p[I_calc],
+            "x"=>x_for_yaml,
+            "p"=>p_for_yaml,
         )
 
         YAML.write_file(path_solution, dict_calc_solution)
diff --git a/test/testcases/pem/MultiVariate_2Normal_3.yml b/test/testcases/pem/MultiVariate_2Normal_3.yml
new file mode 100644
index 0000000..9742565
--- /dev/null
+++ b/test/testcases/pem/MultiVariate_2Normal_3.yml
@@ -0,0 +1,18 @@
+x:
+  -
+    - -1.732050807568877
+    - 0.0
+    - 1.732050807568877
+  -
+    - -1.732050807568877
+    - 0.0
+    - 1.732050807568877
+p:
+  -
+    - 0.16666666666666666
+    - 0.16666666666666666
+    - 0.16666666666666666
+  -
+    - 0.16666666666666666
+    - 0.16666666666666666
+    - 0.16666666666666666
diff --git a/test/testcases/pem/MultiVariate_3Normal_2.yml b/test/testcases/pem/MultiVariate_3Normal_2.yml
new file mode 100644
index 0000000..3f44aba
--- /dev/null
+++ b/test/testcases/pem/MultiVariate_3Normal_2.yml
@@ -0,0 +1,20 @@
+x:
+  -
+    - -1.7320508075688776
+    - 1.7320508075688772
+  -
+    - -1.7320508075688776
+    - 1.7320508075688772
+  -
+    - -1.7320508075688776
+    - 1.7320508075688772
+p:
+  -
+    - 0.16666666666666666
+    - 0.16666666666666669
+  -
+    - 0.16666666666666666
+    - 0.16666666666666669
+  -
+    - 0.16666666666666666
+    - 0.16666666666666669
diff --git a/test/testcases/pem/MultiVariate_3Normal_3.yml b/test/testcases/pem/MultiVariate_3Normal_3.yml
new file mode 100644
index 0000000..fa74aa8
--- /dev/null
+++ b/test/testcases/pem/MultiVariate_3Normal_3.yml
@@ -0,0 +1,26 @@
+x:
+  -
+    - -1.7320508075688776
+    - 0.0
+    - 1.7320508075688772
+  -
+    - -1.7320508075688776
+    - 0.0
+    - 1.7320508075688772
+  -
+    - -1.7320508075688776
+    - 0.0
+    - 1.7320508075688772
+p:
+  -
+    - 0.1666666666666666
+    - -0.0
+    - 0.16666666666666663
+  -
+    - 0.1666666666666666
+    - -0.0
+    - 0.16666666666666663
+  -
+    - 0.1666666666666666
+    - -0.0
+    - 0.16666666666666663
diff --git a/test/testcases/pem/MultiVariate_Normal_2.yml b/test/testcases/pem/MultiVariate_Normal_2.yml
new file mode 100644
index 0000000..766c37a
--- /dev/null
+++ b/test/testcases/pem/MultiVariate_Normal_2.yml
@@ -0,0 +1,8 @@
+x:
+  -
+    - -1
+    - +1
+p:
+  -
+    - 0.5
+    - 0.5