Please suggest the needful.

\[Tau] = 9.2; \[Mu] = 3.5; \[Alpha] = 0.5; \[Gamma] = 0.5; m = 5;
\[Omega] = 1; \[Tau]c = 3; \[Epsilon] = 0.1; \[Eta] = 0.4;
\[CapitalOmega] = 0.1;
a = 4/(9 \[Gamma] \[Mu] ) (\[Mu] Cos[\[Tau] + \[Omega] \[Epsilon]
\[Tau]] + \[Mu] \[Alpha] - \[Eta]);

Tmax = 300;
ClearAll[f, ψ, t]

f[ψ_, τc_] := 2 Ω + (m/
          2)*(-(m/a)*Sin[τc + ω ϵ τc]*
           Sin[ψ]*
           Sin[(τc + ω ϵ τc) - ψ] - (m/a)*
           Sin[τc + ω ϵ τc]*Sin[ψ]*
           Sin[(τc + ω ϵ τc) + ψ]);
tauList = Range[0, 10, 0.05];
psiInit = 1;
forwardData = 
  Reap[Do[sol = 
      NDSolve[{ψ'[t] == f[ψ[t], τ], ψ[0] == 
         psiInit}, ψ, {t, 0, Tmax}, MaxSteps -> Infinity];
     (*remove transient*)
     psiFinal = 
      Mean[Table[ψ[t] /. sol[[1]], {t, Tmax - 50, Tmax, 1}]];
     Sow[{τ, psiFinal}];
     psiInit = psiFinal;  (*continuation*), {τ, tauList}]][[2, 
   1]];
tauListBack = Reverse[tauList];
psiInit = 7;
backwardData = 
  Reap[Do[sol = 
      NDSolve[{ψ'[t] == f[ψ[t], τ], ψ[0] == 
         psiInit}, ψ, {t, 0, Tmax}, MaxSteps -> Infinity];
     psiFinal = 
      Mean[Table[ψ[t] /. sol[[1]], {t, Tmax - 50, Tmax, 1}]];
     Sow[{τ, psiFinal}];
     psiInit = psiFinal;  (*continuation*), {τ, tauListBack}]][[
   2, 1]];
Show[ListLinePlot[forwardData, PlotStyle -> Red], 
 ListLinePlot[backwardData, PlotStyle -> Black], 
 AxesLabel -> {"τc", "ψ"}, PlotRange -> All]