Point charges plot

Here is a beautiful example of using tubes to draw scalar fields (taken from [stackoverflow](https://mathematica.stackexchange.com/questions/687/id-like-to-display-field-lines-for-a-point-charge-in-3-dimensions))
fieldSolve[field_, varlist_, xi0_, tmax_, debug_: False] := Module[
  {xiVec, equationSet, t},
  If[Length[varlist] != Length[xi0], 
   Print["Number of variables must equal number of initial conditions\
\nUSAGE:\n" <> fieldSolve::usage]; Abort[]];
  xiVec = Through[varlist[t]];
  (* Below, Simplify[equationSet] would cost extra time 
    and doesn't help with the numerical solution, so   don't try to simplify. *)

  equationSet = Join[
    Thread[
     Map[D[#, t] &, xiVec] == 
      Normalize[field /. Thread[varlist -> xiVec]]
     ],
    Thread[
     (xiVec /. t -> 0) == xi0
     ]
    ];
  If[debug, 
   Print[Row[{"Numerically solving the system of equations\n\n", 
      TraditionalForm[(Simplify[equationSet] /. t -> "t") // 
        TableForm]}]]];
  (* This is where the differential equation is solved. 
  The Quiet[] command suppresses warning messages because numerical precision isn't crucial for our plotting purposes: *)

  Map[Head, First[xiVec /.
     Quiet[NDSolve[
       equationSet,
       xiVec,
       {t, 0, tmax}
       ]]], 2]
  ]   

fieldLinePlot[field_, varList_, seedList_, opts : OptionsPattern[]] :=
   Module[{sols, localVars, var, localField, plotOptions, 
    tubeFunction, tubePlotStyle, postProcess = {}}, 
   plotOptions = FilterRules[{opts}, Options[ParametricPlot3D]];
   tubeFunction = OptionValue["TubeFunction"];
   If[tubeFunction =!= None,
    tubePlotStyle = Cases[OptionValue[PlotStyle], Except[_Tube]];
    plotOptions = 
     FilterRules[plotOptions, 
      Except[{PlotStyle, ColorFunction, ColorFunctionScaling}]];
    postProcess = 
     Line[x_] :> 
      Join[tubePlotStyle, {CapForm["Butt"], 
        Tube[x, tubeFunction @@@ x]}]
    ];
   If[Length[seedList[[1, 1]]] != Length[varList], 
    Print["Number of variables must equal number of initial \
conditions\nUSAGE:\n" <> fieldLinePlot::usage]; Abort[]];
   localVars = Array[var, Length[varList]];
   localField = 
    ReleaseHold[
     Hold[field] /. 
      Thread[Map[HoldPattern, Unevaluated[varList]] -> localVars]];
   (*Assume that each element of seedList specifies a point AND the \
length of the field line:*)Show[
    Table[
     ParametricPlot3D[
         Evaluate[
          Through[#[t]]], {t, #[[1, 1, 1, 1]], #[[1, 1, 1, 2]]}, 
         Evaluate@Apply[Sequence, plotOptions]
         ] &[fieldSolve[
        localField, localVars, seedList[[i, 1]], seedList[[i, 2]]
        ]
       ] /. postProcess, {i, Length[seedList]}
     ]
    ]
   ];

Options[fieldLinePlot] = 
  Append[Options[ParametricPlot3D], "TubeFunction" -> None];

SetAttributes[fieldSolve, HoldAll];

f2[x_, y_,z_] := {x, y, z}/Norm[{x, y, z}]^3 - ({x, y, z} - {1, 1, 1})/
   Norm[{x, y, z} - {1, 1, 1}]^3

seedList = 
  With[{vertices = .1 N[PolyhedronData["Icosahedron"][[1, 1]]]}, 
   Join[Map[{#, 2} &, vertices], 
    Map[{# + {1, 1, 1}, -2} &, vertices]]];

fieldLinePlot[f2[x, y, z], {x, y, z}, seedList, 
 PlotStyle -> {Orange, Directive["Roughness"->0., "Specularity"->White, "SpecularityIntensity"->10000]}, PlotRange -> All, Lighting->None, Epilog->{
   PointLight[White, 5 {1,1,1}, 100],
   SpotLight[Red, 5 {1,1,1}, 100],
   SpotLight[Green, 5 {1,-1,1}, 100],
   SpotLight[Blue, 5 {-1,1,1}, 100],
   HemisphereLight[LightBlue, Black, 0.1]
  },
 "TubeFunction" -> 
  Function[{x, y, z}, Quiet[Clip[Norm[f2[x, y, z]], {2, 40}]/300]], 
 Background -> Black, ImageSize->500]