Basically I am generating a grid in space (so a vector where v[[1]]=position 1 in x) and a random number of particles in space (so a vector where p[[1]]= position of particle 1 on space). Afterwards, I want to see in which spaces of the grid the particles are and do somethings. The Do is the specific part of the code I want to optimize. rho = Compile[{{xp, _Real, 1}, {xg, _Real, 1}, {carga, _Real}, {np, _Real}, {ng, _Real}, {dx, _Real}}, n = ng + 1; ρ = Table[carga*np/(dx*ng), n]; ρ[[ng + 1]] = 0; Do[If[Abs[xp[[i]] - xg[[j]]] <= dx, ρ[[j]] = ρ[[j]] - carga*(1 - (Abs[xp[[i]] - xg[[j]]])/dx)/dx; r += 1, 0] , {i, 1, np}, {j, 1, ng + 1}]; ρ[[1]] = ρ[[1]] + ρ[[ng + 1]];] Basically I am worried that I am thinking too much in terms of other languages and maybe mathematica has a better way to do this instead of using If[Abs[xp[[i]] - xg[[j]]] <= dx That manually makes a lot of operations