Okay, I went ahead and looked at performance a bit. First, here's an alternate to your rho function:

rho2[xp_, xg_, carga_, np_, ng_, dx_] :=
 Module[
  {rhoTable = Table[carga*np/(dx*ng), 1 + ng]},
  rhoTable[[1 + ng]] = 0;
  With[
   {neighbors = Nearest[xg -> {"Index", "Distance"}, xp, {All, dx}]},
   Map[(rhoTable[[#[[1]]]] -= (carga*(1 - (#[[2]])/dx)/dx)) &, 
    neighbors, {-2}];
   rhoTable[[1]] += rhoTable[[ng + 1]];
   {Length[Flatten@neighbors]/2, rhoTable}]]

I used local variables rather than update global ones. It looked like you were interested in number of updates, so I returned both that and the rho table, although I don't know if number of updates would ever be anything other than 2*nparticles. Here's the timing I got on my machine for your version:

(result1 = rho[x, xgrid, q, nparticles, ngrid, dx]; {r, \[Rho]}) // AbsoluteTiming // Short
(*{0.11498, {4000, {-12152.9, 8155.76, -4859.24, <<195>>, -11347.9, 468.88, -60255.7}}}*)

And here's the timing with the version using Nearest:

(result2 = rho2[x, xgrid, q, nparticles, ngrid, dx]) // AbsoluteTiming // Short
(*{0.02478, {4000, {-12152.9, 8155.76, -4859.24, <<195>>, -11347.9, 468.88, -60255.7}}}*)

And a quick check:

result1[[2]] == result2[[2]]
(*True*)

Diogo,

I am confused. so you are saying my code should be something like this?

Almost. However, Do[] loops are generally not recommended in Mathematica. Also, the FractionalPart is always positive so the Abs is not necessary. I would use Map or MapThread instead of the loop.

For example, your posted code takes 0.2 seconds to run compiled. This code gets the same result but takes only 0.007 seconds on my M1 machine without compilation:

AbsoluteTiming[rr = Table[q*nparticles/(dx*ngrid), n];
 rr[[ngrid + 1]] = 0;
 MapThread[(rr[[#1]] = rr[[#1]] - q*(1 - (#2))/dx) &, {pos, frac}];
 MapThread[(rr[[#1]] = rr[[#1]] - q*((#2))/dx) &, {pos + 1, frac}];
 rr[[1]] = rr[[1]] + rr[[ngrid + 1]];]

Regards,

Neil

Diogo,

What you seem to be doing is finding where each point lies on a grid. There is a much simpler way to do this rather than to compare the value to every point on the grid. I would first compute the integer part of the x values to find where on the grid you start:

pos = Floor[x/dx] + 1;

I did not use Ceil because I think you always round up even if the error is zero (but if I'm wrong you can change this)

Next get the fractional part:

frac = FractionalPart[x/dx];

or get both at the same time with MixedFractionParts[].

With this information you can set your rho values with

 MapThread[rhoFormula,{pos,frac}]

in which you index into rho using pos and do your equivalent computation of the Do loop using frac.This works because you do all this looping to find the position in the grid just below and just above your x values. pos is the position in the grid below your value and pos+1 is the position in the grid above your value. frac and (1 - frac) gives you the expression Abs[xp[[i]] - xg[[j]]])/dx which is used in your calculations.

I hope this helps.

Regards,

Neil

Just a clarification... Is it expected that each position of \[Rho] can get updated multiple times during each inner loop?