Solving a system of non-linear equations and obtaining a phase portrait

Still a newbie in Mathematica, I used streamplot function to generate a bunch of streamlines for the non-linear system below. How do I solve this system of equations numerically and generate the appropriate solution curves?

y'[t]=m[x,y,c,k]  

x'[t]=n[x,y,c]

Where

m[x_, y_, c_, k_] := (((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]) + k
n[x_, y_, c_] := -((2 (x - 1) (y - 1))/((y - 1)^2 + x^2))*((1 - Exp[-(((y - 1)^2 + x^2)/c)])/(((y - 1)^2 + x^2)/c)).

c ranges from 0.01 to 100, and k ranges from 0 to 1.

Thank you in advance!

Answer

eq1 = x'[t] == -((2 (x[t] - 1) (y[t] - 1))/((y[t] - 1)^2 + x[t]^2))*((1 - 
        Exp[-(((y[t] - 1)^2 + x[t]^2)/c)])/(((y[t] - 1)^2 + x[t]^2)/c));
eq2 = y'[t] == (((y[t] - 1)^2 - x[t]^2)/((y[t] - 1)^2 + x[t]^2)*((1 - 
          Exp[-(((y[t] - 1)^2 + x[t]^2)/c)])/(((y[t] - 1)^2 + x[t]^2)/
          c))) - (ExpIntegralE[1, ((y[t] - 1)^2 + x[t]^2)/c]) + k;

c = 3; k = 0.5;

sol[x0_?NumericQ] := First@NDSolve[{eq1, eq2, x[0] == x0, y[0] == x0}, {x, y}, {t, 0, 200}]


pp = ParametricPlot[Evaluate[{x[t], y[t]} /. sol[#] & /@ Range[0.3, 1, 0.2]], {t, 0,200}];

sp = StreamPlot[{-((2 (x - 1) (y - 1))/((y - 1)^2 + x^2))*((1 - 
         Exp[-(((y - 1)^2 + x^2)/c)])/(((y - 1)^2 + x^2)/
         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]) + k}, {x, -5, 1}, {y, 0, 2}];

Show[pp, sp]

sol1 = ParametricNDSolveValue[{eq1, eq2, x[0] == x0, y[0] == y0}, {x, y}, {t, -100, 100}, 
{x0, y0}];

points = Join[Table[{-2, y}, {y, -5, 5, 0.1}], Table[{2, y}, {y, -5, 5, 0.1}]];  

ParametricPlot[sol1 @@@ points//Evaluate, {t, -100, 100}, PlotRange -> {{-5, 5}, {-5, 5}}]