diff --git a/src/CartesianProduct.jl b/src/CartesianProduct.jl
index 5be6984be9..ac6fd53ea4 100644
--- a/src/CartesianProduct.jl
+++ b/src/CartesianProduct.jl
@@ -6,7 +6,8 @@ export CartesianProduct,
        array,
        block_indices,
        variable_indices,
-       block_structure
+       block_structure,
+       same_block_structure
 
 """
     CartesianProduct{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}
diff --git a/src/concrete_intersection.jl b/src/concrete_intersection.jl
index 23305f54d8..0da0b6f59e 100644
--- a/src/concrete_intersection.jl
+++ b/src/concrete_intersection.jl
@@ -780,7 +780,7 @@ function intersection(rm::ResetMap{N, <:AbstractPolytope},
 end
 
 """
-function intersection(X::CartesianProductArray{N},
+        intersection(X::CartesianProductArray{N},
                       Y::AbstractPolyhedron{N}) where {N}
 
 Return the intersection of the cartesian product of a finite number of convex sets and a polyhedron.
@@ -800,26 +800,29 @@ This function takes into account that if the polyhedron is unbounded, the inters
 is only needed to be taken in the elements of the cartesian product array (subsets of
 variables, or "blocks") which are constrained; those which are not constrained do not need to be intersected.
 """
-function intersection_combine(X::CartesianProductArray{N},
+function intersection(X::CartesianProductArray{N},
                       Y::AbstractPolyhedron{N},
                       blocks::Dict{Int,Int},
                       oa=Hyperrectangle) where {N}
 
     low_set, vars, block_structure = initialize_vars(X, blocks)
 
-    low_intersection = intersection(approx_low_set, Approximations.project(Y,vars, LinearMap))
+    low_intersection = intersection(low_set, Approximations.project(Y, vars, LinearMap))
 
     if isempty(low_intersection)
         return EmptySet()
     end
 
+    println(block_structure)
     decomposed_low_set = Approximations.decompose(low_intersection, block_structure, oa)
 
-    return combine_resulted_cpa(decomposed_low_set, X)
+    return combine_resulted_cpa(decomposed_low_set, X, blocks)
 end
 
-function combine_resulted_cpa(l_cpa::CartesianProductArray{N}, o_cpa::CartesianProductArray{N}) where {N}
-    result = CartesianProductArray(length(X.array), N)
+function combine_resulted_cpa(l_cpa::CartesianProductArray{N},
+                              o_cpa::CartesianProductArray{N},
+                              blocks::Dict{Int,Int}) where {N}
+    result = CartesianProductArray(length(l_cpa.array), N)
 
     index = 1
     for bi in 1:length(o_cpa.array)
@@ -842,7 +845,7 @@ function initialize_vars(X::CartesianProductArray{N}, blocks::Dict{Int,Int}) whe
         n = dim(X.array[bi])
         push!(low_set.array, X.array[bi])
         append!(vars, variable_indices(X, bi, blocks[bi]))
-        push!(block_structure, start_index:(start_index + n))
+        push!(block_structure, start_index:(start_index + n - 1))
         start_index += n
     end
 
@@ -854,7 +857,7 @@ end
 
 function intersection(X::CartesianProductArray{N},
                       Y::AbstractPolyhedron{N},
-                      oa::Union{Type{<:Approximations.AbstractDirections},Type{<:Hyperrectangle}, Type{<:Interval}}=Hyperrectangle) where {N}
+                      oa=Hyperrectangle) where {N}
 
     if isbounded(Y)
         # no free variables
@@ -869,7 +872,7 @@ end
 
 # symmetric method
 function intersection(Y::AbstractPolyhedron{N}, X::CartesianProductArray{N},
-                      oa::Union{Type{<:Approximations.AbstractDirections},Type{<:Hyperrectangle}, Type{<:Interval}}=Hyperrectangle) where {N}
+                      oa=Hyperrectangle) where {N}
     intersection(X, Y, oa)
 end
 
@@ -881,3 +884,33 @@ end
 function intersection(X::CartesianProductArray{N}, U::Universe{N}) where {N<:Real}
     return intersection(U, X)
 end
+
+"""
+        intersection(X::CartesianProductArray{N}, Y::CartesianProductArray{N})
+
+Return the intersection between cartesian products of a finite number of convex sets.
+
+### Input
+
+ - `X` -- cartesian product of a finite number of convex sets
+ - `Y` -- cartesian product of a finite number of convex sets
+
+### Output
+
+The decomposed set which represents concrete intersection between `X` and `Y`
+
+### Algorithm
+
+This algorithm intersect corresponding blocks between sets.
+"""
+function intersection(X::CartesianProductArray{N}, Y::CartesianProductArray{N}) where {N<:Real}
+    @assert same_block_structure(array(X), array(Y)) "block structure has to be the same"
+
+    result = CartesianProductArray(length(X.array), N)
+
+    for i in 1:length(X.array)
+        push!(result, intersection(X.array[i], Y.array[i]))
+    end
+
+    return result
+end