If the derivative odefun(t,y) depends on the four "control" functions of t (or even of both t and y), then these functions would be incorporated into the definition of yprime = odefun(t,y). I don't see how the fact that the solver is fixed-step or variable-step should matter, except that fixed-step simulations are numerically suspect to begin with. Even small amounts of instability will be magnified when integrating for any length using the (Forward) Euler method with fixed step sizes.
ODE23 provides a dizzying array of options, including event functions, nonnegativity contraints, and output functions. You can even have it return a struct (use one output) and then pass that struct to ODEXTEND to extend the integration. Using that facility, you could conceivably implement a pseudo-fixed-size integration with a loop, but the only benefit I can think of to that is if the control functions change in a non-smooth fashion at discrete points. In that case it would be best not to try to integrate over these points as if they weren't there, rather to step on, and subsequently from, them precisely.