Error : no method matching getindex

Your problem is println("X[$i] = ", value(X[i])). The capital X is a single variable, not a vector.

I’d refactored your code a little, but you still have a problem: it doesn’t converge because z isn’t part of the optimization model (and so it can pick any value).

using JuMP
import Juniper
import Ipopt

function inner_loop(old_x, old_y, R)
    m = 15 # num customers
    NumVar = 5 # num hospitals
    a = [0, 2.5, -2.5, 2, -2.5] # x cordinates of the circle
    b = [0, 0.3, 0.25, -2.5, -2.25]  # y cordinates of the circle
    xs = [0, 2.5, -2.5, 2, -2.5] # x cordinates of the circle
    ys = [0, 0.3, 0.25, -2.5, -2.25]  # y cordinates of the circle
    r = [2, 2, 2, 2, 2] # Radiuus of circle
    R = 0.0
    coverage = [
        1 0 0 0 1;
        0 0 0 1 1;
        1 0 0 0 1;
        0 0 1 1 1;
        0 0 0 1 1;
        0 1 0 0 1;
        0 0 1 1 0;
        0 0 0 1 0;
        1 0 0 0 0;
        1 1 0 0 0;
        0 1 0 0 0;
        0 0 1 0 0;
        0 0 1 0 0;
        1 0 0 0 0;
        0 0 0 1 0
    ]
    for i in 1:length(a)
        R = max(R, sqrt((a[i] - old_x)^2 + (b[i] - old_y)^2) + r[i])
    end    
    model = Model(
        optimizer_with_attributes(
            Juniper.Optimizer,
            "nl_solver" => optimizer_with_attributes(
                Ipopt.Optimizer, 
                MOI.Silent() => true,
            ),
        ),
    )
    @variable(model, x[1:NumVar], Bin)
    @variable(model, R >= 0)
    @variable(model, X >= 0)
    @variable(model, Y >= 0)
    @variable(model, z >= 0)
    @objective(model, Min, R)
    @constraint(model, [i=1:m], coverage[i, :]' * x >= 1)
    @NLconstraint(
        model, 
        [i=1:NumVar],
        (sqrt((xs[i] - X)^2 + (ys[i] - Y)^2) + r[i]) * x[i] <= R
    )
    optimize!(model)
    for i in 1:NumVar
        println("x[$i] = ", value(x[i]))
    end
    return value(X), value(Y), value(z), value(R)
end

function main()
    x, y, z, R = 0.0, 0.0, Inf, 0.0
    while abs(x^2 + y^2 - z) >= 0.01
        x, y, z, R = inner_loop(x, y, R)
    end
    print("The abs value is: abs(x^2+y^2-z) = " , abs(x^2 + y^2 - z))
    return x, y, z, R
end

main()