ClearAll[ifnPart];
ifnPart["Domain"] = Sequence[1];   (*bounding box for domain*)
ifnPart["X1"] = Sequence[1, 1];    (* lower bound for first coordinate *)
ifnPart["X2"] = Sequence[1, 2];    (* upper bound for first coordinate *)
ifnPart["Version"] = Sequence[2, 1]; 
ifnPart["Flags"] = Sequence[2, 2]; (*flags indicating properties:
    bit field positions - inferred, perhaps mistaken
      $extrapolation=0; whether to warn about extrapolation
      $fullArrayBit=1; interpolation data is a full array (not ragged)
      $packed=2; packed array form (???)
      $repeatedBit=4; whether repeated abscissae are permitted*)
ifnPart["DerivativeOrder"] = Sequence[2, 3]; (*max derivative order*)
ifnPart["NGrid"] = Sequence[2, 4];  (*number of points in each coordinate grid*)
ifnPart["InterpolationOrder"] = Sequence[2, 5]; (*interpolation order*)
ifnPart["Derivative"] = Sequence[2, 6];  (*derivative to evaluate:0-->f[x], 1-->f'[x],...*)
ifnPart["Periodic"] = Sequence[2, 7];
(*ifnPart["??"]=Sequence[2,8];*)
(*ifnPart["??"]=Sequence[2,9];*)
ifnPart["ExtrapolationHandler"] = Sequence[2, 10];
(*ifnPart["??"]=Sequence[2,11];*)
(*ifnPart["??"]=Sequence[2,12];*)
(*ifnPart["??"]=Sequence[2,13];*)
ifnPart["Coordinates"] = Sequence[3];  (*list of lists, abscissae of interpolation grid*)
ifnPart["InterpolationData"] = Sequence[4]; (*interpolation data (values or coefficients)*)
ifnPart["Offsets"] = Sequence[4, 2]; (*offsets in function/derivative array (PackedArrayForm)*)
ifnPart["FlatData"] = Sequence[4, 3]; (*flattened function/derivative values (PackedArrayForm)*)
ifnPart["InterpolationStructure"] = Sequence[5]; (*{Automatic}, or dense output interpolation structure:
   list of types for each unit/subinterval*)
ifnPart["UnitIndices"] = Sequence[5, 1, 1]; (*dense output:
   Indices (to grid) for corresponding coefficients*)
ifnPart["UnitTypes"] = Sequence[5, 1, 2];   (*dense output types:
   Automatic | NDSolve`CubicHermite | NDSolve`LocalSeries | ChebyshevT*)

ifnPart["Properties"] = 
 Cases[DownValues[ifnPart], Verbatim[ifnPart][prop_] :> prop, Infinity];

ifnPart["ValidPartQ", "Chebyshev" | "Local Taylor Series", "UnitIndices" | "UnitTypes", _] := True;
ifnPart["ValidPartQ", _, "UnitIndices" | "UnitTypes", _] := False;
ifnPart["ValidPartQ", "Hermite", "Offsets" | "FlatData", Developer`PackedArrayForm] := True;
ifnPart["ValidPartQ", _, "Offsets" | "FlatData", _] := False;
ifnPart["ValidPartQ", method_String, part_String, _] /; 
 MemberQ[method, "Chebyshev" | "Local Taylor Series" | "Hermite"] &&
  MemberQ[part, ifnPart["Properties"]] := True;
ifnPart["ValidPartQ", _, _, _] := False;
ifnPart[if_InterpolatingFunction, part_String] /; 
   ifnPart["ValidPartQ", if["InterpolationMethod"], part, if[[4, 1]]] :=
  if~Part~ifnPart[part];

ClearAll[ifnPart];
ifnPart["Domain"] = Sequence[1];   (*bounding box for domain*)
ifnPart["X1"] = Sequence[1, 1];    (* lower bound for first coordinate *)
ifnPart["X2"] = Sequence[1, 2];    (* upper bound for first coordinate *)
ifnPart["Version"] = Sequence[2, 1]; 
ifnPart["Flags"] = Sequence[2, 2]; (* flags indicating properties:
    bit field positions - inferred, perhaps mistaken
      $extrapolation=0; whether to warn about extrapolation
      $fullArrayBit=1; interpolation data is a full array (not ragged)
      $machPrecBit=2; whether the data (f,f',...) is MachinePrecision
      $repeatedBit=4; whether repeated abscissae are permitted *)
ifnPart["DerivativeOrder"] = Sequence[2, 3]; (*max derivative order*)
ifnPart["NGrid"] = Sequence[2, 4];  (*number of points in each coordinate grid*)
ifnPart["InterpolationOrder"] = Sequence[2, 5]; (*interpolation order*)
ifnPart["Derivative"] = Sequence[2, 6];  (*derivative to evaluate:0-->f[x], 1-->f'[x],...*)
ifnPart["Periodic"] = Sequence[2, 7];
(*ifnPart["??"]=Sequence[2,8];*)
(*ifnPart["??"]=Sequence[2,9];*)
ifnPart["ExtrapolationHandler"] = Sequence[2, 10];
(*ifnPart["??"]=Sequence[2,11];*)
(*ifnPart["??"]=Sequence[2,12];*)
(*ifnPart["??"]=Sequence[2,13];*)
ifnPart["Coordinates"] = Sequence[3];  (*list of lists, abscissae of interpolation grid*)
ifnPart["InterpolationData"] = Sequence[4]; (*interpolation data (values or coefficients)*)
ifnPart["Offsets"] = Sequence[4, 2]; (*offsets in function/derivative array (PackedArrayForm)*)
ifnPart["FlatData"] = Sequence[4, 3]; (*flattened function/derivative values (PackedArrayForm)*)
ifnPart["InterpolationStructure"] = Sequence[5]; (*{Automatic}, or dense output interpolation structure:
   list of types for each unit/subinterval*)
ifnPart["UnitIndices"] = Sequence[5, 1, 1]; (*dense output:
   Indices (to grid) for corresponding coefficients*)
ifnPart["UnitTypes"] = Sequence[5, 1, 2];   (*dense output types:
   Automatic | NDSolve`CubicHermite | NDSolve`LocalSeries | ChebyshevT*)

ifnPart["Properties"] = 
 Cases[DownValues[ifnPart], Verbatim[ifnPart][prop_] :> prop, Infinity];

ifnPart["ValidPartQ", "Chebyshev" | "Local Taylor Series", "UnitIndices" | "UnitTypes", _] := True;
ifnPart["ValidPartQ", _, "UnitIndices" | "UnitTypes", _] := False;
ifnPart["ValidPartQ", "Hermite", "Offsets" | "FlatData", Developer`PackedArrayForm] := True;
ifnPart["ValidPartQ", _, "Offsets" | "FlatData", _] := False;
ifnPart["ValidPartQ", method_String, part_String, _] /; 
 MemberQ[method, "Chebyshev" | "Local Taylor Series" | "Hermite"] &&
  MemberQ[part, ifnPart["Properties"]] := True;
ifnPart["ValidPartQ", _, _, _] := False;
ifnPart[if_InterpolatingFunction, part_String] /; 
   ifnPart["ValidPartQ", if["InterpolationMethod"], part, if[[4, 1]]] :=
  if~Part~ifnPart[part];

This handles the "Local Taylor series" and "Chebyshev" interpolation methods, but not the packed or unpacked "Hermite" ones. One thing I would recommend, and implemented below, is not to worry about cutting off the interpolation exactly at tmin or tmax. Rather, take the smallest interval containing tmin and tmax with endpoints on the interpolation grid. Actually the code below might include one more grid point at each end than is necessary. This adds very little overhead. On the other hand, dealing with dividing a subinterval with these methods might be a much bigger pain than it's worth. In any case, the function below preserves these (more accurate!) interpolation methods. It seems a shame to throw away the extra accuracy. One might do it if speed is more important than accuracy.

ClearAll[ifTake];
dupeLast[list_] := Append[list, Last@list];
ifTake[if_InterpolatingFunction, {tmin_?NumericQ, tmax_?NumericQ}] /; 
    MatchQ[if["InterpolationMethod"], "Local Taylor series" | "Chebyshev"] :=
  Module[{coords, newif = Hold @@ if, span},
   coords = First@if["Coordinates"];
   span = Clip[
     SparseArray[UnitStep[coords - tmin] UnitStep[tmax - coords]][
        "AdjacencyLists"][[{1, -1}]] + {-1, 1}, {1, Length@coords}];
   newif[[1]] = {coords[[span]]};
   (*newif[[2,{2,3}]]={??};*) (* might need to update flags/deriv. order? *)
   newif[[2, 4]] = 1 + Differences@span;
   newif[[3]] = Developer`ToPackedArray@{coords[[Span @@ span]]};
   newif[[4]] = Join[
     {if[[4, First@span, 1 ;; 2]]},
     if[[4, First@span + 1 ;; Last@span]]
     ];
   newif[[5, 1, 1]] = Join[
     {{1}},
     1 + if[[5, 1, 1, First@span + 1 ;; Last@span]] - 
       if[[5, 1, 1, First@span, -1]] // dupeLast
     ];
   newif[[5, 1, 2]] = Join[
     {Automatic},
     if[[5, 1, 2, First@span + 1 ;; Last@span]] // dupeLast
     ];
   InterpolatingFunction @@ newif
   ];

OP's first example:

if1 = n /. 
   NDSolve[{n'[t] == n[t] (1 - n[t - 2]), n[0] == 0.1}, 
     n, {t, 0, 20000}, MaxSteps -> ∞][[1]];
if1["InterpolationMethod"]
(*  "Local Taylor series"  *)

ifTake[if1, {0.5, 1.5}] // AbsoluteTiming
(*  {0.028446, InterpolatingFunction[{{0.478428, 1.53819}}, <>]}  *)

ifTake[if1, {1000.5, 1010.5}] // AbsoluteTiming
(*  {0.03142, InterpolatingFunction[{{1000.46, 1010.55}}, <>]}  *)

Chebyshev example:

if3 = NDSolveValue[{y''[x] + y[x] == 0, y[0] == 1, y'[0] == 0}, 
   y, {x, 0, 100}, InterpolationOrder -> All, Method -> "Extrapolation"];
if3["InterpolationMethod"]
(*  "Chebyshev"  *)

ifTake[if3, {10.5, 11.5}]
(*  InterpolatingFunction[{{9.03213, 12.7461}}, <>]  *)

This handles the "Local Taylor series", "Chebyshev" and both the packed and unpacked "Hermite" interpolation method. One thing I would recommend, and implemented below, is not to worry about cutting off the interpolation exactly at tmin or tmax. Rather, take the smallest interval containing tmin and tmax with endpoints on the interpolation grid. Actually the code below might include one more grid point at each end than is necessary. This adds very little overhead. On the other hand, dealing with dividing a subinterval with these methods might be a much bigger pain than it's worth. In any case, the function in the "Code dump" below preserves these (more accurate!) interpolation methods. It seems a shame to throw away the extra accuracy. One might do it if speed is more important than accuracy.

OP's first example:

if1 = n /. NDSolve[{n'[t] == n[t] (1 - n[t - 2]), n[0] == 0.1}, 
     n, {t, 0, 20000}, MaxSteps -> ∞][[1]];
if1["InterpolationMethod"]
(*  "Local Taylor series"  *)

ifTake[if1, {0.5, 1.5}] // AbsoluteTiming
(*  {0.028446, InterpolatingFunction[{{0.478428, 1.53819}}, <>]}  *)

ifTake[if1, {1000.5, 1010.5}] // AbsoluteTiming
(*  {0.03142, InterpolatingFunction[{{1000.46, 1010.55}}, <>]}  *)

Chebyshev example:

if3 = NDSolveValue[{y''[x] + y[x] == 0, y[0] == 1, y'[0] == 0}, 
   y, {x, 0, 100}, InterpolationOrder -> All, Method -> "Extrapolation"];
if3["InterpolationMethod"]
(*  "Chebyshev"  *)

ifTake[if3, {10.5, 11.5}]
(*  InterpolatingFunction[{{9.03213, 12.7461}}, <>]  *)

Code dump

InterpolatingFunction parts. Execute

ifnPart["Properties"] (* OR *)
DownValues[ifnPart][[;; Length@ifnPart["Properties"] + 1]]

to see a list of the part names. Some are valid only for the "Hermite" method, some for "Chebyshev" or "Local Taylor Series". The function of some of the parts are as yet unknown to me. It as complete as my current knowledge, and more than is needed. But it seemed worth sharing.

ClearAll[ifnPart];
ifnPart["Domain"] = Sequence[1];   (*bounding box for domain*)
ifnPart["X1"] = Sequence[1, 1];    (* lower bound for first coordinate *)
ifnPart["X2"] = Sequence[1, 2];    (* upper bound for first coordinate *)
ifnPart["Version"] = Sequence[2, 1]; 
ifnPart["Flags"] = Sequence[2, 2]; (*flags indicating properties:
    bit field positions - inferred, perhaps mistaken
      $extrapolation=0; whether to warn about extrapolation
      $fullArrayBit=1; interpolation data is a full array (not ragged)
      $packed=2; packed array form (???)
      $repeatedBit=4; whether repeated abscissae are permitted*)
ifnPart["DerivativeOrder"] = Sequence[2, 3]; (*max derivative order*)
ifnPart["NGrid"] = Sequence[2, 4];  (*number of points in each coordinate grid*)
ifnPart["InterpolationOrder"] = Sequence[2, 5]; (*interpolation order*)
ifnPart["Derivative"] = Sequence[2, 6];  (*derivative to evaluate:0-->f[x], 1-->f'[x],...*)
ifnPart["Periodic"] = Sequence[2, 7];
(*ifnPart["??"]=Sequence[2,8];*)
(*ifnPart["??"]=Sequence[2,9];*)
ifnPart["ExtrapolationHandler"] = Sequence[2, 10];
(*ifnPart["??"]=Sequence[2,11];*)
(*ifnPart["??"]=Sequence[2,12];*)
(*ifnPart["??"]=Sequence[2,13];*)
ifnPart["Coordinates"] = Sequence[3];  (*list of lists, abscissae of interpolation grid*)
ifnPart["InterpolationData"] = Sequence[4]; (*interpolation data (values or coefficients)*)
ifnPart["Offsets"] = Sequence[4, 2]; (*offsets in function/derivative array (PackedArrayForm)*)
ifnPart["FlatData"] = Sequence[4, 3]; (*flattened function/derivative values (PackedArrayForm)*)
ifnPart["InterpolationStructure"] = Sequence[5]; (*{Automatic}, or dense output interpolation structure:
   list of types for each unit/subinterval*)
ifnPart["UnitIndices"] = Sequence[5, 1, 1]; (*dense output:
   Indices (to grid) for corresponding coefficients*)
ifnPart["UnitTypes"] = Sequence[5, 1, 2];   (*dense output types:
   Automatic | NDSolve`CubicHermite | NDSolve`LocalSeries | ChebyshevT*)

ifnPart["Properties"] = 
 Cases[DownValues[ifnPart], Verbatim[ifnPart][prop_] :> prop, Infinity];

ifnPart["ValidPartQ", "Chebyshev" | "Local Taylor Series", "UnitIndices" | "UnitTypes", _] := True;
ifnPart["ValidPartQ", _, "UnitIndices" | "UnitTypes", _] := False;
ifnPart["ValidPartQ", "Hermite", "Offsets" | "FlatData", Developer`PackedArrayForm] := True;
ifnPart["ValidPartQ", _, "Offsets" | "FlatData", _] := False;
ifnPart["ValidPartQ", method_String, part_String, _] /; 
 MemberQ[method, "Chebyshev" | "Local Taylor Series" | "Hermite"] &&
  MemberQ[part, ifnPart["Properties"]] := True;
ifnPart["ValidPartQ", _, _, _] := False;
ifnPart[if_InterpolatingFunction, part_String] /; 
   ifnPart["ValidPartQ", if["InterpolationMethod"], part, if[[4, 1]]] :=
  if~Part~ifnPart[part];

Taking part of an InterpolatingFunction:

ClearAll[ifTake];
dupeLast[list_] := Append[list, Last@list];
iDataTake["Local Taylor series" | "Chebyshev", data_, span_] := Join[
   {data[[First@span, 1 ;; 2]]}, data[[First@span + 1 ;; Last@span]]
   ];
iDataTake["Hermite", data : {Developer`PackedArrayForm, _, _}, 
   span : {s1_, s2_}] := ReplacePart[
   data,
   {Rest@{ifnPart["Offsets"]} ->
     data[[2, s1 ;; s2 + 1]] - data[[2, s1]],
    Rest@{ifnPart["FlatData"]} ->
     data[[3, data[[2, s1]] + 1 ;; data[[2, s2 + 1]] ]]}
   ];
iDataTake["Hermite", data : {__List}, span_] := data[[Span @@ span]];
iStructureTake["Local Taylor series" | "Chebyshev", structure_, 
   span_] := ReplacePart[structure,
   {Rest@{ifnPart["UnitIndices"]} ->
     Join[
      {{1}},
      1 + structure[[##2 &@ifnPart["UnitIndices"], First@span + 1 ;; Last@span]] - 
        structure[[##2 &@ifnPart["UnitIndices"], First@span, -1]] // 
       dupeLast
      ], 
    Rest@{ifnPart["UnitTypes"]} ->
     Join[
      {Automatic},
      structure[[##2 &@ifnPart["UnitTypes"], 
         First@span + 1 ;; Last@span]] // dupeLast
      ]}
   ];
iStructureTake["Hermite", structure_, span_] := structure;
ifTake[if_InterpolatingFunction, {tmin_?NumericQ, tmax_?NumericQ}] /; 
   Length@if["Domain"] == 1 :=
  Module[{coords, newif = Hold @@ if, span, method},
   method = if["InterpolationMethod"];
   coords = First@if["Coordinates"];
   span = Clip[
     SparseArray[UnitStep[coords - tmin] UnitStep[tmax - coords]][
        "AdjacencyLists"][[{1, -1}]] + {-1, 1}, {1, Length@coords}];
   
   newif[[ifnPart["Domain"]]] =
    {coords[[span]]};
   (*newif[[ifnPart[if,"Flags"]]]=??;
     newif[[ifnPart[if, "DerivativeOrder"]}]]=??;*) (* not needed?? *)
   newif[[ifnPart["NGrid"]]] =
    1 + Differences@span;
   newif[[ifnPart["Coordinates"]]] =
    Developer`ToPackedArray@{coords[[Span @@ span]]};
   newif[[ifnPart["InterpolationData"]]] =
    iDataTake[method, if[[ifnPart["InterpolationData"]]], span];
   newif[[ifnPart["InterpolationStructure"]]] =
    iStructureTake[method, if[[ifnPart["InterpolationStructure"]]], span];
    
   InterpolatingFunction @@ newif
   ];