Its going to be something like this: (if I understand what you want)

\[CapitalGamma] = 5/3;
c4 = 6;
smax = 2;
Rinit = 1.05;
\[Gamma]1 = 0.8;
\[Gamma]2 = 0.81;
\[Gamma]uinit = (\[Gamma]1 + \[Gamma]2)/2;
sols = Table[
  Reap@NDSolve[
    smax = 2; {D1[s] == 
      2 (-1 + R[s]) R[
        s] (-(-1 + \[CapitalGamma]) \[Gamma]u[s] + 
         c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
             2) (-2 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
            s]^(2 - \[CapitalGamma]) + 
         c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
             2) (-1 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
            s]^-\[CapitalGamma]), 
     N1[s] == -(1 + \[Gamma]u[s]^2) (1 - \[CapitalGamma] + 
         c4 (-1 + R[s])^(1/2 - \[CapitalGamma]/2) R[
            s]^(-(3/2) (-1 + \[CapitalGamma])) (2 + 
            4 R[s] (-1 + \[CapitalGamma]) - 
            3 \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
            s]^(1 - \[CapitalGamma])), 
     R'[s] == D1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
     Derivative[1][\[Gamma]u][s] == N1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
     R[0] == Rinit, \[Gamma]u[0] == \[Gamma]uinit, 
     WhenEvent[{D1[s] == 0, N1[s] == 0}, 
      Sow[{s, D1[s], N1[s], Abs[N1[s]/D1[s]]}]; smax = s;
      If[Abs[N1[s]/D1[s]] > 1, {\[Gamma]uinit = \[Gamma]2, 
        Print[\[Gamma]uinit]}];
      If[Abs[N1[s]/D1[s]] < 1, {\[Gamma]uinit = \[Gamma]1, 
        Print[\[Gamma]uinit]}]]}, {R, \[Gamma]u}, {s, 0, smax}], 
  5]; 
data =  Map[({R[s], \[Gamma]u[s]} /. #) &, 
      sols[[All, 1, 1]]]; 
ParametricPlot[data, {s, 0, smax}, 
     AspectRatio -> 1, AxesLabel -> {"R", "\[Gamma]u"}]

The idea is to make a list of solutions and then plot the list of results.

Regards,

Neil

Your bug was that you included the line smax = s in the WhenEvent[]. By doing this you forced the integration to stop just when you wanted the If statement to alter your constants.

This works fine:

\[CapitalGamma] = 5/3;
c4 = 6;
smax = 2;
Rinit = 1.05;
\[Gamma]1 = 0.8;
\[Gamma]2 = 0.81;
\[Gamma]uinit = (\[Gamma]1 + \[Gamma]2)/2;
{sol1, values} = 
  Reap@NDSolve[{D1[s] == 
      2 (-1 + R[s]) R[
        s] (-(-1 + \[CapitalGamma]) \[Gamma]u[s] + 
         c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
             2) (-2 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
            s]^(2 - \[CapitalGamma]) + 
         c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
             2) (-1 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
            s]^-\[CapitalGamma]), 
     N1[s] == -(1 + \[Gamma]u[s]^2) (1 - \[CapitalGamma] + 
         c4 (-1 + R[s])^(1/2 - \[CapitalGamma]/2) R[
            s]^(-(3/2) (-1 + \[CapitalGamma])) (2 + 
            4 R[s] (-1 + \[CapitalGamma]) - 
            3 \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
            s]^(1 - \[CapitalGamma])), 
     R'[s] == D1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
     Derivative[1][\[Gamma]u][s] == N1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
     R[0] == Rinit, \[Gamma]u[0] == \[Gamma]uinit, 
     WhenEvent[{D1[s] == 0, 
       N1[s] == 0}, {Sow[{s, D1[s], N1[s], Abs[N1[s]/D1[s]]}],
       If[Abs[N1[s]/D1[s]] > 1, \[Gamma]uinit = \[Gamma]2; 
        Print[\[Gamma]uinit]];
       If[Abs[N1[s]/D1[s]] < 1, \[Gamma]uinit = \[Gamma]1; 
        Print[\[Gamma]uinit]]}]}, {R, \[Gamma]u}, {s, 0, smax}];

Regards,

Neil

Agapi,

There is no need to loop. WhenEvent has the ability to change states when events happen and continue the integration. Look at the example in the documentation of WhenEvent. The first example is exactly what you need. I did not dive into your code to see what state you are checking in the If statement but I'll assume its [Gamma]u. You would do something like

WhenEvent[N1[s]* D1[s] == 0,If[\[Gamma]u > 1, \[Gamma]2 = \[Gamma]uinit,\[Gamma]1 = \[Gamma]uinit]]

Although I really do not understand what you are doing so you will have to change the action above to fit your situation. The documentation examples are good and you can execute multiple statements inside a list in the WhenEvent.

Regards

Neil