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 Here is a part of the code I am using eps = 8.85*10^(-12); q = 80; m = 50; L = 2; nparticles = 2000; ngrid = 200; vth = 2; u = 1; dx = L/ngrid; x = RandomVariate[UniformDistribution[{0, L}], nparticles]; r = 0; xgrid = Range[0, L, dx] // N; rho[x, xgrid, q, nparticles, ngrid, dx]; // AbsoluteTiming r