The following Julia program uses Makie to create the image shown below, a "cube with holes." However, I would like to eliminate the holes, and also make the cube nearer the shape of an actual cube, with only slightly rounded corners and edges. I've tried changing the exponent to values other than ^2 and also passing the range directly to the volume! function, to no avail. Examples of the types of cubes desired can be found in three-rounded-box images at CodeSandbox.

using Makie, GLMakie
fig = Figure()
range = LinRange(-1, 1, 100) # 100-element LinRange{Float64, Int64}
cube = [ (x.^2 + y.^2 + z.^2) for x = range, y = range, z = range ] # 100×100×100 Array{Float64, 3}
ax = Axis3( fig[1,1], aspect = :data, azimuth = 1.17 * pi, viewmode = :fit, title = "Cube" )
volume!( cube, algorithm = :iso, isorange = 0.05, isovalue = 1.7 )
fig

Like this?

using Makie, GLMakie
fig = Figure()
range = LinRange(-1, 1, 100) # 100-element LinRange{Float64, Int64}
cube = [ (abs.(x).^10 + abs.(y).^10 + abs.(z).^10) for x = range, y = range, z = range ] # 100×100×100 Array{Float64, 3}
ax = Axis3( fig[1,1], aspect = :data, azimuth = 1.17 * pi, viewmode = :fit, title = "Cube" )
volume!( cube, algorithm = :iso, isorange = 0.05, isovalue = 1 )
fig

Or you can use a parametric surface, a supershape:

using Makie, GLMakie

function r(phi; a, b, m, n1, n2, n3)
    return 1 / (abs(cos(m*phi/4)/a)^n2 + (abs(sin(m*phi/4)/b)^n3))^(1/n1) 
end

phi = (-pi/2):0.01:(pi/2)
theta = (-pi):0.01:(pi)
x = [
    r(theta; a=1, b=1, m=4, n1=10, n2=10, n3=10) * cos(theta) *
    r(phi;   a=1, b=1, m=4, n1=10, n2=10, n3=10) * cos(phi) 
    for phi in phi, theta in theta
]
y = [
    r(theta; a=1, b=1, m=4, n1=10, n2=10, n3=10) * sin(theta) *
    r(phi;   a=1, b=1, m=4, n1=10, n2=10, n3=10) * cos(phi) 
    for phi in phi, theta in theta
]
z = [
    r(phi;   a=1, b=1, m=4, n1=10, n2=10, n3=10) * sin(phi) 
    for phi in phi, theta in theta
]
fig, _ = surface(x, y, z)
fig

Better to use a function:

function r(phi; a, b, m, n1, n2, n3)
    return 1 / (abs(cos(m*phi/4)/a)^n2 + (abs(sin(m*phi/4)/b)^n3))^(1/n1) 
end

function supershape(p1, p2)
    phi = (-pi/2):0.01:(pi/2)
    theta = (-pi):0.01:(pi)
    x = [
        r(theta; a=p1.a, b=p1.b, m=p1.m, n1=p1.n1, n2=p1.n2, n3=p1.n3) * cos(theta) *
        r(phi;   a=p2.a, b=p2.b, m=p2.m, n1=p2.n1, n2=p2.n2, n3=p2.n3) * cos(phi) 
        for phi in phi, theta in theta
    ]
    y = [
        r(theta; a=p1.a, b=p1.b, m=p1.m, n1=p1.n1, n2=p1.n2, n3=p1.n3) * sin(theta) *
        r(phi;   a=p2.a, b=p2.b, m=p2.m, n1=p2.n1, n2=p2.n2, n3=p2.n3) * cos(phi) 
        for phi in phi, theta in theta
    ]
    z = [
        r(phi;   a=p2.a, b=p2.b, m=p2.m, n1=p2.n1, n2=p2.n2, n3=p2.n3) * sin(phi) 
        for phi in phi, theta in theta
    ]
    return (x = x, y = y, z = z)
end

params1 = (a=1, b=1, m=4, n1=10, n2=10, n3=10)
params2 = params1
x, y, z = supershape(params1, params2)

fig, _ = surface(x, y, z)

EDIT

The rendering is not nice with surface. It's better to do a mesh:

function parametricMesh(f, umin, umax, vmin, vmax, nu, nv)
    u_ = LinRange(umin, umax, nu)
    v_ = LinRange(vmin, vmax, nv)
    vertices = Array{Float64}(undef, nu*nv, 3)
    triangles = Array{Int64}(undef, 2*(nu-1)*(nv-1), 3)
    k = 1
    for i in 1:nv
        v = v_[i]
        for j in 1:nu
            vertices[k,:] = f(u_[j], v)
            k = k+1
        end
    end
    k = 1
    for i in 1:(nv-1)
        for j in 1:(nu-1)
            a = (i-1) * nu + j
            b = (i-1) * nu + j + 1
            c = i*nu + j + 1
            d = i*nu + j
            triangles[2*(k-1)+1,:] = [b, a, d]
            triangles[2*k,:] = [c, b, d]
            k = k+1
        end
    end
    return (vertices = vertices, triangles = triangles)
end

function r(phi; a, b, m, n1, n2, n3)
    return 1 / (abs(cos(m*phi/4)/a)^n2 + (abs(sin(m*phi/4)/b)^n3))^(1/n1) 
end

p1 = (a=1, b=1, m=4, n1=10, n2=10, n3=10)
p2 = p1

function f(phi, theta)
    x = r(theta; a=p1.a, b=p1.b, m=p1.m, n1=p1.n1, n2=p1.n2, n3=p1.n3) * cos(theta) *
        r(phi;   a=p2.a, b=p2.b, m=p2.m, n1=p2.n1, n2=p2.n2, n3=p2.n3) * cos(phi)
    y = r(theta; a=p1.a, b=p1.b, m=p1.m, n1=p1.n1, n2=p1.n2, n3=p1.n3) * sin(theta) *
        r(phi;   a=p2.a, b=p2.b, m=p2.m, n1=p2.n1, n2=p2.n2, n3=p2.n3) * cos(phi)
    z = r(phi;   a=p2.a, b=p2.b, m=p2.m, n1=p2.n1, n2=p2.n2, n3=p2.n3) * sin(phi)
    return [x, y, z]
end

vertices, triangles = parametricMesh(f, -pi/2, pi/2, -pi, pi, 50, 50)

mesh(vertices, triangles, color = "yellow")