I don't know the cause of these jumps when finding a solution with a discontinuous reference function

I have a simple function with error dynamics edot = -b0*e, which means that if I track a reference function Fref, the output F has state equation Fdot = -Frefdot - b0*(Fref - F). This can be solved using a discrete approximation, or using a two-state augmentation to get rid of the input derivative in the state equation. The problem occurs when Fref is discontinuous, e.g. if Fref = square(t).

MATLAB's ODE solvers are not intended to work with discontinuous functions, as previously pointed out in this question/answer. However, I'm working with a controller that provides discontinuous inputs, so I need an ok solution (i.e. not high fidelity, but reliably consistent) solution to this problem.

I tried to write a short, very simple 4th order Runge-Kutta solver so that I could solve the system above. The problem is that it still is inconsistent, and since it's a square wave and my RK4 is just arithmetic within the function at this point, that makes no sense to me (the behavior at time t = .5 should match that at time t = 1.5). Can someone please explain why my solver isn't consistent over the length of the square wave?

I attached my files (a script and two short functions). Running the solver_Analysis.m script shows results. The figure below shows the output. The reference is plotted, as is a discrete solution (to show what I expect), ode45 (to show why I can't use that; ode23 works slightly better but still with errors), and my RK4 solver.