hor[x_, y_, c_] := (((y - 1)^2 - x^2)/((y - 1)^2 + x^2)*((
  1 - Exp[-(((y - 1)^2 + x^2)/c)])/(((y - 1)^2 + x^2)/
  c))) - (ExpIntegralE[1, ((y - 1)^2 + x^2)/c]);
ver[x_, y_, c_] := -(((y - 1)^2 - x^2)/((y - 1)^2 + x^2)*((
   1 - Exp[-(((y - 1)^2 + x^2)/c)])/(((y - 1)^2 + x^2)/
   c))) - (ExpIntegralE[1, ((y - 1)^2 + x^2)/c]);
mix[x_, y_, c_] := -((2 x (y - 1))/((y - 1)^2 + x^2))*((
1 - Exp[-(((y - 1)^2 + x^2)/c)])/(((y - 1)^2 + x^2)/c));

dim[x_, y_, c_, w_] := 1/2 (ver[x, y, c] - hor[x, y, c]) + w;

top := dim[x, y, c, w] - Sqrt[((dim[x, y, c, w])^2 + (mix[x, y, c])^2)]
bottom := mix[x, y, c]
points = Join[Table[{-3, y}, {y, -5, 5, 0.5}], Table[{3, y}, {y, -5, 5, 0.5}]];
param = {w -> 0.5, c -> 1};
odesystem = {x'[t], y'[t]} == ({bottom, top} /. {p : x | y :> p[t]} /. param);

paramsol = ParametricNDSolveValue[{odesystem, {x[0], y[0]} == {m, k}}, {x, y}, {t, -100, 100}, {m, k},
Method -> {StiffnessSwitching,
 Method -> {ExplicitRungeKutta, Automatic}}, PrecisionGoal -> 10, AccuracyGoal -> 8, MaxSteps -> 100000, MaxStepSize -> 0.1];
soln = paramsol @@@ points;
Graphics[{Black, Arrowheads[0], Arrow[Transpose@Through[#["ValuesOnGrid"]] & /@ soln],{Black, Thick, Line[{{0, 0.5}, {0, 1.5}}]}}, PlotRange -> 5, Frame ->    True,PlotRangePadding -> Scaled[.02]]