I am trying to solve friedman acceleration equation and continuity equation numerically using Runge Kutta Dormand Prince Method . My code was doing excellent for solving coupled differential equations. But, when I am trying to solve this equation , I am getting wrong results , I am getting same variation of scale factor for both matter dominated and radiation dominated case while one should vary as a^{2/3} and other should vary as a^{1/2}. Can someone please tell me the possible reasons why this is happen ?

\[\Ddota= -1/6(\rho + 3p)\]