One option is to use NDSolve components.
Also interpolation itself is not that time consuming (you can set InterpolationOrder -> 1).
Values on grid can be then extracted from InterpolatingFunction with DifferentialEquations`InterpolatingFunctionAnatomy` package.
Take a look at tutorial/NDSolveStateData and tutorial/NDSolvePackages.