I really don't understand what you are trying to do. Can you give a simple example? At first I thought you wanted to hit an event and then change some settings and continue. Next you said that you only wanted to do one simulation and stop at the event, reset everything, change a setting and start over. This is what your code is currently doing. Four of the simulations are identical because once you make your setting change the answer is always the same. Your code runs in both cases and produces results so without really understanding what you want, it's hard to be more helpful.

Regards

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