Is there a way to sort NSolve solution (roots) automatically?

OK, I've had this problem before myself. Here's a solution to a harder problem? Take your code to generate the roots list.

Remove["Global`*"];//Quiet
ClearAll;
blackbox[a_]=(x-Sin[a]-1)*(x-Sin[a-Pi]-1)
(-1+x-Sin[a]) (-1+x+Sin[a])
xx=Range[0,2*Pi,0.01];
roots=Table[x/.NSolve[blackbox[b]==0,x],{b,xx}];

Add a complication by just scrambling the results.

rootScramble = Map[#[[RandomSample[Range[2]]]] &, roots];

Plot a subset. Ugly!

Now, making the assumption that the results are smooth in the zeroth and 1st derivative, look at each point pair, and then permute the next pair in line to minimize a Norm involving both their values and derivatives. The idea is that while in the "point" space they might have the same value, in the phase space of both point and derivative, they differ and can be discriminated.

First create the permutator, which takes in as the first argument (assuming two curves, but should work for any number) {{p1,p2},{del1,del2}}. The dels are finite differences between the that point and the last one, stored up as the algorithm runs. The second argument is the point that needs permuting, {q1,q2}.

permMatch[{vals_, dels_}, l2_] := 
 Module[{l1 = {vals, dels}, ps, bestPerm},

  ps = Map[({#, # - vals}) &, Permutations[l2]];

  bestPerm = Sort[ps, Norm[l1 - #1] < Norm[l1 - #2] &] // First
  ]

The output of permMatch is a permutation {{qx,qy},{qx-p1,qy-p2}} that is "closest" to {{p1,p2},{del1,del2}}.

Now just use FoldList to process rootScramble, initializing with the {{0,0},{0,0}} to get it started then dropping it at the end.

res = FoldList[permMatch, {{0, 0}, {0, 0}}, rootScramble]//Rest;

Discard the differences and transpose

finalResult = Map[(# // First) &, res] // Transpose;

Plot the full before and after.

rootScramble // Transpose // ListLinePlot

finalResult // ListLinePlot

This could easily be extended to include the 2nd derivative as a discriminator, so for example where two curves touch and have the same value and slope, it could still discriminate them.

EDIT

Applied it your problem above...

e5 = {e1, e2, e3, e4};
roots = Transpose[{e1, e2, e3, e4}];
ListPlot[e5, PlotRange -> {-0.2, 0.2}]

res = FoldList[permMatch, {{0, 0, 0, 0}, {0, 0, 0, 0}}, roots] // Rest;
finalResult = Map[(# // First) &, res] // Transpose;
finalResult // ListLinePlot

EDIT (last one, I promise)

Tried this on a much messier problem

f1 = Sin[t] + 3 Cos[t] + 3 Sin[2 t]
f2 = 2 Sin[t] + 3 Cos[4 t] - 3 Sin[2 t]
f3 = 2 Sin[3 t] + 3 Cos[4 t] - 3 Sin[5 t] + 1
f4 = 2 Sin[6 t] + 3 Cos[5 t] - 3 Sin[t]

Now scramble it all up.

    scrambled = Table[{f1, f2, f3, f4}[[RandomSample[Range[4]]]], {t, 0, 2π, .05}]

When I ran the algorithm...

Not a right match. Needed to weight the derivative portion more, and now it works out.

WEIGHT = 10;    

permMatchWeight[{vals_, dels_}, l2_] := 
 Module[{l1 = {vals, dels}, ps, bestPerm}, 
  ps = Map[({#, WEIGHT (# - vals)}) &, Permutations[l2]];
  bestPerm = Sort[ps, Norm[l1 - #1] < Norm[l1 - #2] &] // First]

A standard way to solve an equation $f(x,y)=0$ for the function $y(x)$ satisfying a condition $y(x_0)=y_0$ is to integrate the differentiated equation as an IVP: $$f_x(x,y) + f_y(x,y)\,y' = 0,\quad y(x_0) = y_0$$ This will work on equations with multiple solutions as long as the solution functions do not cross or touch each other. If there are such singular points, then one may "blow up" the singularity by differentiating the equation further and using the derivatives of $y(x)$ to separate the solutions. Simple crossings may be resolved by differentiating once. The resulting differential equation is second-order, and the dependent variables $\langle y(x),\,y'(x)\rangle$ will be distinct for distinct solutions. For solutions that are tangent, higher-order derivatives would be needed depending on the order of tangency.

eqn = (y - Cos[x] + 2) (y + Cos[x] + 2) (y - Cos[x] - 2) (2 y + Sin[x] y - 2) == 0;
eqnx = eqn /. y -> y[x];  (* equation with y as a function of x *)
ode1 = D[eqnx, x];        (* first-order ODE *)
ode2 = D[ode1, x];        (* second-order ODE *)

ics0 = NSolve[eqnx /. x -> -2 Pi, y[-2 Pi]];
ics1 = First@Solve[ode1, y'[x]] /. x -> -2 Pi /. ics0;
ics = Apply[Equal, Join[ics0, ics1, 2], {2}]
(*
{{y[-2 \[Pi]] == -3., Derivative[1][y][-2 \[Pi]] == 0.},
 {y[-2 \[Pi]] == -1., Derivative[1][y][-2 \[Pi]] == 0.},
 {y[-2 \[Pi]] == 1., Derivative[1][y][-2 \[Pi]] == -0.5},
 {y[-2 \[Pi]] == 3., Derivative[1][y][-2 \[Pi]] == 0.}}
*)

NDSolveValue[{ode2, #}, y, {x, -2 Pi, 2 Pi}, 
    Method -> {"FixedStep", 
      Method -> {"Projection", "Invariants" -> ({eqnx, ode1} /. Equal -> Subtract)}}, 
    StartingStepSize -> 1/100] & /@ ics // ListLinePlot

Remark on robustness: One should note that I used a "FixedStep" method, together with the "Projection" method. The projection method projects the solution step back onto the original equations, basically using Newton's method. This ensures the accuracy of the steps. However, Newton's method tends not to work well near the singular points, that is, the crossings. The fixed-step method was used to step over the crossings and avoid the numerical difficulties they present. With adaptive step-size methods, the step size tends to diminish near the crossings and trouble is encountered. Also, the wrong StartingStepSize, such as 2 Pi/100 in this example, leads to a step landing on a singularity and fails. Some adjustment may be necessary in some problems.