How to find a numerical antiderivative with NIntegrate methods?

We construct an NDSolve method which can pass an NIntegrate method to NIntegrate to set up an integration rule. We define a method nintegrate implements such a method. The requirements are

• the ODE is of the form y'[x] == f[x], and
• the NIntegrate method returns an interpolatory rule.

Example:

foo = NDSolveValue[{y'[x] == Sin[x^2], y[0] == 0}, y, {x, 0, 15}, 
  Method -> nintegrate, InterpolationOrder -> All]

Error plot:

Plot[
 Evaluate@RealExponent[Integrate[Sin[x^2], x] - foo[x]],
 {x, 0, 15},
 GridLines -> {Flatten@foo@"Grid", None}, (* show steps *)
 PlotRange -> {-18.5, 0.5}]

Another example:

foo = NDSolveValue[{y'[x] == Sin[x^2], y[0] == 0}, y, {x, 0, 15}, 
  Method -> {nintegrate, 
    Method -> {"ClenshawCurtisRule", "Points" -> 33}}, 
  InterpolationOrder -> All, WorkingPrecision -> 32, 
  PrecisionGoal -> 24, MaxStepFraction -> 1, StartingStepSize -> 15]

Error plot:

Block[{$MaxExtraPrecision = 500},
 ListLinePlot[
  Integrate[Sin[x^2], x] - foo[x] /. x -> Subdivide[0, 15, 1000] // 
   RealExponent, DataRange -> {0, 15}, PlotRange -> {-35.5, 0.5}, 
  GridLines -> {Flatten@foo@"Grid", None}]
 ]

Code for method

nintegrate::nintode = 
  "Method nintegrate requires an ode of the form ``'[``] == f[``]";
nintegrate::nintinit = 
  "NIntegrate method `` did not return an interpolatory integration rule.";

nintegrate[___]["StepInput"] = {"F"["T"], "H", "T", "X", "XP"};
nintegrate[___]["StepOutput"] = {"H", "XI"};
nintegrate[rule_, order_, ___]["DifferenceOrder"] := order;
nintegrate[___]["StepMode"] := Automatic

Options@nintegrate = {Method -> "GaussKronrodRule"};
getorder[points_, method_] :=
  Switch[method
   , "GaussKronrodRule" | "GaussKronrod",
   (* check points should be odd ??? *)
   With[{gp = (points - 1)/2},
    If[OddQ[gp], 3 gp + 2, 3 gp + 1]
    ]
   , "LobattoKronrodRule",
   (* check points should be odd ??? *)
   With[{glp = (points + 1)/2},
    If[OddQ[glp], 3 glp - 2, 3 glp - 3]
    ]
   , "GauseBerntsenEspelidRule",
   2 points - 1
   , "NewtonCotesRule",
   If[OddQ[points], points, points - 1]
   , _, points - 1
   ];
nintegrate /: 
  NDSolve`InitializeMethod[nintegrate, stepmode_, sd_, rhs_, state_, 
   mopts : OptionsPattern[nintegrate]] := 
  Module[{prec, order, norm, rule, xvars, tvar, imeth},
   xvars = NDSolve`SolutionDataComponent[state@"Variables", "X"];
   tvar = NDSolve`SolutionDataComponent[state@"Variables", "T"];
   If[Length@xvars != 1,
    Message[nintegrate::nintode, First@xvars, tvar, tvar];
    Return[$Failed]];
   If[! VectorQ[rhs["FunctionExpression"][
       N@NDSolve`SolutionDataComponent[sd, "T"],
       Sequence @@ xvars],
      NumericQ
      ],
    Message[nintegrate::nintode, First@xvars, tvar, tvar];
    Return[$Failed]];
   prec = state@"WorkingPrecision";
   norm = state@"Norm";
   
   imeth = Replace[Method /. mopts, Automatic -> "GaussKronrodRule"];
   rule = 
    NIntegrate[1, {x, 0, 1}, 
     Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0, 
       Method -> imeth},
     WorkingPrecision -> prec,
     IntegrationMonitor :>
      (Return[Through[#@"GetRule"], NIntegrate] &)];
   rule = Replace[rule, {
      {(NIntegrate`GeneralRule | NIntegrate`ClenshawCurtisRule)[idata_]} :>
       idata,
      _NIntegrate :>
       Return[$Failed],
      _ :>  (* What happened here? *)
       (Message[nintegrate::nintinit, Method -> imeth];
        Return[$Failed])
      }];
   order = 
    getorder[Length@First@rule, imeth /. {m_String, ___} :> m];
   
   nintegrate[rule, order, norm]
   ];

(rule : nintegrate[int_, order_, norm_, ___])[
   "Step"[rhs_, h_, t_, x_, xp_]] := 
  Module[{prec, tt, xx, dx, normh, err, hnew, temp},
   (* Norm scaling will be based on current solution y. *)
   normh = (Abs[h] temp[#1, x] &) /. {temp -> norm};
   tt = Rescale[int[[1]], {0, 1}, {t, t + h}];
   xx = rhs /@ tt;
   dx = h*int[[2]].xx;
   (* Compute scaled error estimate *)
   err = h*int[[3]].xx // normh;
   hnew = Which[
     err > 1 (* Rejected step: reduce h by half *)
     , dx = $Failed; h/2
     , err < 2^-(order + 2), 2 h
     , err < 1/2, h
     , True, h Max[1/2, Min[9/10, (1/(2 err))^(1/(order + 1))]]
     ];
   (* Return step data along with updated method data *)
   {hnew, dx}];