Actually no I am not sure, and as I am new in Mathematica, I really don't know if I have any other way to solve it! I would appreciate it if you could give me any hint.

Your code has some syntax problems. This makes sense:

\[Epsilon]0 = 8.85*10^-12; \[Epsilon]g = 
 42*\[Epsilon]0; \[Xi] = 0.1; R1 = 50*10^-9;
\[Rho]ag = 10.49*10^3; cd = 0.5 10^-8;

\[Epsilon]rel = ((\[Epsilon]g - \[Epsilon]0)/(\[Epsilon]g + \
\[Epsilon]0));

x1 = 500 10^-9; z1 = 5 10^-9; z1s = -5 10^-9; zs = 
 5 10^-9; x2 = -500 10^-9; z2 = 5 10^-9; z2s = -5 10^-9;


U[xs_] := (1/(Sqrt[(xs - x1)^2 + (zs - z1)^2])^3 - (3*(xs - 
           x1)^2)/(Sqrt[(xs - x1)^2 + (zs - z1)^2])^5 + 
     1/(Sqrt[(xs - x2)^2 + (zs - z2)^2])^3 - (3*(xs - 
           x2)^2)/(Sqrt[(xs - x2)^2 + (zs - 
             z2)^2])^5) - \[Epsilon]rel*(1/(Sqrt[(xs - x1)^2 + (zs - 
              z1s)^2])^3 - (3*(xs - 
            x1)^2)/(Sqrt[(xs - x1)^2 + (zs - z1s)^2])^5 + 
      1/(Sqrt[(xs - x2)^2 + (zs - z2s)^2])^3 - (3*(xs - 
            x2)^2)/(Sqrt[(xs - x2)^2 + (zs - z2s)^2])^5 + 
      1/(Sqrt[(zs - zss)^2])^3);

Felectrical = -U'[xs[t]];
Fsurfaceviscosity = -cd*xs'[t];
ma = 4/3*\[Pi]*\[Xi]*(R1)^3*\[Rho]ag*xs''[t];

Ftotal = Felectrical + Fsurfaceviscosity - ma;

NDSolve[{Ftotal == 0, xs[0] == 100 10^-9, 
  xs[2700] == 300 10^-9}, xs, {t, 0, 3031}]

but it may not be what you meant. However, it does not compute, because of numerical problems.