In this paper, a newly proposed embedded Runge–Kutta fourth order with four stages arithmetic Heronian mean method is employed and verified in determining the distribution of thermodynamic variables inside the extra-solar protoplanets formed through gravitational instability at their initial stages. In specific, the case of conduction–radiation is considered regarding the transference of heat inside the protoplanets. A general brief theoretical framework for the proposed numerical method is stated in addition to pseudo code followed by error estimation description. The results based on newly proposed explicit RKAHeM(4,4) method are found to be optimal and efficient in comparison with the ones obtained with classical fourth order Runge–Kutta method.
