How can I vectorize this code in R? Maybe with the apply() function?
First of all it should be noted that the code you posted does not actually replicate the output of the dist function, because the line:

distancematrix[i, j] &lt;- sum(abs(myMatrix[i,] - myMatrix[j,]))


does not calculate the Euclidean distance; it should be:

distancematrix[i, j] &lt;- sqrt(sum((myMatrix[i,] - myMatrix[j,]) ^ 2))


Here are two solutions that rely on apply. They are simplified, and in particular do not take advantage of the symmetry of the distance matrix (which, if considered, would lead to a 2-fold speedup). First, generate some test data:

# Number of data points
N &lt;- 2000
# Dimensionality
d &lt;- 10
# Generate data
myMatrix = matrix(rnorm(N * d), nrow = N)


For convenience, define:

# Wrapper for the distance function
d_fun &lt;- function(x_1, x_2) sqrt(sum((x_1 - x_2) ^ 2))


The first approach is a combination of apply and sapply:

system.time(
    D_1 &lt;- 
        apply(myMatrix, 1, function(x_i) 
            sapply(1:nrow(myMatrix), function(j) d_fun(x_i, myMatrix[j, ]))
    )
)

   user  system elapsed 
 14.041   0.100  14.001 


while the second uses only apply (but going over the indices, which are paired using expand.grid):

system.time(
    D_2 &lt;- 
        matrix(apply(expand.grid(i = 1:nrow(myMatrix), j = 1:nrow(myMatrix)), 1, function(I) 
            d_fun(myMatrix[I[["i"]], ], myMatrix[I[["j"]], ])
        )
    )
)

   user  system elapsed 
 39.313   0.498  39.561 


However, as expected both are much slower than dist:

system.time(
    distancematrix &lt;- as.matrix(
        dist(myMatrix, method = "euclidean", diag = T, upper = T)
    )
)

   user  system elapsed 
  0.337   0.054   0.388