Performing multi-input GP regression using Stheno.jl

I am attempting to perform multi-input GP regression using Stheno.jl but unable to find any example and keep getting dimension errors.

What I want to learn is a multidimensional function y = g(x,z) using GP.

I attempted to convert the single-dimensional examples as:

using AbstractGPs
using LinearAlgebra
using Stheno
using Plots
using Random

# Pre-defining the hyper-parameters
l1 = 0.4
s1 = 0.2
# l2 = 5.0
# s2 = 1.0

# Defining GPPP(GaussianProcessProbabilisticProgramme) object 
# f = @gppp let
#     f1 = s1 * stretch(GP(Matern52Kernel()), 1 / l1)
#     f2 = s2 * stretch(GP(SEKernel()), 1 / l2)
#     f3 = f1 + f2
# end;

 # Defining GPPP(GaussianProcessProbabilisticProgramme) object  with θ as hyperpramter vector
function build_model(θ::NamedTuple)
    return @gppp let
        f1 = θ.s1 * stretch(GP(SEKernel()), 1 / θ.l1) # stretch here is used to push the length scale in ZeroMean GP
#         f2 = θ.s2 * stretch(GP(SEKernel()), 1 / θ.l2)
#         f3 = f1 + f2
    end
end


function g(x,z)  
           x.^2 +2x + x.^3 +z
end

 Random.seed!(3)
points = randn(10,1)
 Random.seed!(2)
z_points = randn(10,1)
# x = GPPPInput(:f, points)
const y = vec(g(points,z_points))


const x = [GPPPInput(:f1, vec(points)) GPPPInput(:f1, vec(z_points))];# this is to convert inputs into the GPPP input types

# σ²_n = 0.02;
# fx = f(x, σ²_n);
# const y = rand(fx);


using ParameterHandling
using ParameterHandling: value, flatten

# Defining a NamedTuple 
θ = (
    # Short length-scale and small variance.
    l1 = positive(0.4),
    s1 = positive(0.2),



    # Observation noise variance -- we'll be learning this as well. Constrained to be
    # at least 1e-3.
    s_noise = positive(0.1, exp, 1e-3),
)
θ_flat_init, unflatten = flatten(θ);# Making a vector out of the NamedTupel
unpack = value ∘ unflatten;

Is there any example doing the same?