What is the interpolation function used to interpolate the ODE solution using deval?

Hi all,

I am currently making a custom ODE solver, and one functionality which I find useful in default ode solvers, and I would like to immitate in my case, is that if one writes their function as

sol = ode45(odefun,tspan,y0)

then sol can be used with deval in order to interpolate the solution at a specific user-defined time vector, after the integration.

With that said I would like to ask, what is the interpolation function used to perform the above task? Is it for example a cubic, or a spline interpolation? Does it depend on the ODE solver used (my custom solver is similar to ode23tb)?

I have looked into the MATLAB ode functions but haven't managed to find an answer in the above.

Thanks for your help in advance.