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 t2/3 and other should vary as t1/2. Can someone please tell me the possible reasons why this is happen ?