I am using sparse array as variables in the following:

using Ipopt,JuMP
model = Model(Ipopt.Optimizer)
n = 4
@variable(model, X[i=1:n, j=1:i-1] >=0)
@variable(model, y[1:n]>=0)
@variable(model, z[1:n])
@NLconstraint(model, sum(y[i]*sum(X[i,j]*z[j] for j in 1:i-1) for i in 1:n) == 1/6 );
optimize!(model)
println()
@show value.(X)
@show value.(y);

@show defaults to showing only 6 digits, but the rest of them are still there. Just use e.g. println(value.(X)) and println(value.(Y)).

(That being said, I’m not sure what the default tolerances are in JuMP and Ipopt, but I doubt it’s converging to anything close to 16 accurate digits by default. Accuracy ≠ precision.)