InterpolatingFunction with growing Domain

I want to find roots of a function that is very slow to calculate (the function itself involves a bunch of FindMaxima) so I interpolate the function in a region where I think it likely the the root is, since I have very little certainty where the root is I have to interpolate an absurdly large region.

Is there a way to have InterpolatingFunction remember the original function and expand its domain when evaluation is attempted outside it?

-
If you are willing to try different interpolation methods then there may be something here of relevance. –  Daniel Lichtblau Apr 4 '13 at 20:09
@DanielLichtblau That's some really neat ways to do interpolation, thanks! –  ssch Apr 4 '13 at 20:23

Here is a way that sadly it has a little bit of overhead making it slower to evaluate. When argument is outside domain it interpolates towards that direction and creates a new InterpolatingFunction by merging together with the old one.

(*   f: function to interpolate
x0: Starting position
size: minimum interval to interpolate each time
order: InterpolatingOrder
*)
AutoInterpolatingFunction[f_, x0_, size_, opts:OptionsPattern[FunctionInterpolation]] :=
Module[{
wrap,
if},
wrap["if"] =
FunctionInterpolation[
f[\[FormalX]], {\[FormalX], x0 - size, x0 + size},
opts];
wrap["f"] = f;
wrap["xmin"] = x0 - size;
wrap["xmax"] = x0 + size;
(* Argument in domain *)
wrap[x_?NumericQ] := wrap["if"][x] /; wrap["xmin"] <= x <= wrap["xmax"];
(* Interpolate a bit further *)
wrap[x_?NumericQ] :=
Module[{newf},
If[x > wrap["xmax"],
newf = FunctionInterpolation[
wrap["f"][\[FormalX]],
{\[FormalX], wrap["xmax"], Max[wrap["xmax"] + size, x]},
opts];
wrap["if"] =
Interpolation[{Join[wrap["if"]["Grid"], Rest@newf["Grid"]],
Join[wrap["if"]["ValuesOnGrid"],
Rest@newf["ValuesOnGrid"]]}\[Transpose]];
wrap["xmax"] = wrap["if"]["Domain"][[1, 2]]
,
newf =
FunctionInterpolation[
wrap["f"][\[FormalX]], {\[FormalX],
Min[x, wrap["xmin"] - size], wrap["xmin"]},
opts];
wrap["if"] =
Interpolation[{Join[newf["Grid"], Rest@wrap["if"]["Grid"]],
Join[newf["ValuesOnGrid"],
Rest@wrap["if"]["ValuesOnGrid"]]}\[Transpose]];
wrap["xmin"] = wrap["if"]["Domain"][[1, 1]]
];
wrap["if"][x]];
wrap
]


A simple speed test:

aip = AutoInterpolatingFunction[Sin[#] Cos[#] &, 0, 1];
ip = FunctionInterpolation[Sin[x] Cos[x], {x, 0, 1}];
AbsoluteTiming[Do[ip[0.5], {100000}]]
AbsoluteTiming[Do[aip[0.5], {100000}]]
(* 1.7s vs 0.7s *)

-