Like most things in science, the answer to "how do I model x" is "it depends how accurate you want to be." The simplest thing would be to pick some number for what you think the surface temperature of the planet should be, and then use the Stefan-Boltzmann law to find out how much the surface radiates. Of course, the problem is that we likely don't know the planet's surface temperature. If you want your model to only be for soon after the planet forms, we can do pretty well by assuming the planet is a sphere of uniform temperature caused by all of the gravitational binding energy.
If it supposed to apply for a time long after the formation of the planet, looking at surface temperatures for other planets in our solar system for inspiration might be a good starting point. However, you have to keep in mind that you'll probably want to look at temps for the dark side of the planet, because trying to solve this for planets where one side is in the sun makes it significantly harder as the symmetry of the problem is destroyed. So, for the rest of the answer I'm going to assume the planet in question is far enough away from its star that there isn't such more light coming from the direction of it's star than there is from any other direction in space.
What follows is a more accurate, but much more complicated model. Be warned: this is a fairly sophisticated mathematical problem and requires mathematical experience up to at least partial differential equations. I'm not sure what your mathematical experience is, so if you expand on that in the comments I can add a section more tailored to you.
With that out of the way, I think L.Dutch has an answer that's on the right track-- this is just attempting to flesh it out a bit more. First off, to model how the heat moves through Earth, we use the heat equation:
$$\frac{\partial u}{\partial t} = \alpha \nabla^2u + \frac{q}{c_p \rho}$$
where $u$ is the heat density, $c_p$ the specific heat capacity, $\rho$ the density, and q is the heat generated\lost in the material due to outside sources. Note that the temperature is related to these quantities by $T=\frac{u}{c_p \rho}$.
Now, technically this is all you need, along with some boundary values. In practice, it might be unclear what some of these values should correspond to, so I'll walk through that. First, it will help simplify stuff if we assume spherical symmetry. This isn't a requirement, and if we want to model a planet close to it's star or be super duper accurate we can't do it. But, since we assumed it wasn't, and since this is already complicated enough, this allows us to make the substitution
$$\nabla^2 \rightarrow \frac{1}{r^2}\frac{\partial}{\partial r}(r^2 \frac{\partial}{\partial r})$$
and the assumption that $u = u(r, t)$ (ie there is no $\theta$ or $\phi$ dependence). Now, to avoid boundary condition problems, we will make the spatial domain of our problem all of $\mathbb{R}^3$. But since this means the problem includes space, which can't conduct heat, we must make sure that $\alpha(r)=0$ when $r$ is greater than the radius of our planet.
Next, we have to determine the form for $q$. This will either make complete sense or none at all and require a ton of explanation, so I'll just leave the result which is that
$$q=-4\pi R_p^2 \sigma \epsilon (\frac{u}{c_p \rho})^4\delta(r-R_p) + \mathbf{1}_{r \in [0, R_p]}D_0r^2e^{-t/\lambda}$$
where the left term accounts for radiative cooling of the planet and the right for radioactive generation of heat within the planet as per L.Dutch's answer. All the variable names are the standard ones from the Stefan-Boltzmann law, $R_p$ is the radius of the planet, and $\mathbf{1}$ is the indicator function.
Finally all that's left is to specify boundary conditions, which should just specify $u(r,0)$ and that $u$ is bounded. The exact form of these will be that the temperature of space (ie $r>R_p$) will be the temperature of the cosmic microwave background, and that inside it will be some distribution of the initial binding energy of the planet (uniform is easiest and probably pretty accurate). I think these conditions should be enough to solve the equation, but it's possible that the non-linear $u^4$ term messes stuff up and the equation is no longer parabolic, in which case I don't have enough math knowledge to help.
After all that, you "simply" have to solve this monstrosity of an equation. A numerical solver is probably easiest, but you could attempt an analytical solution. Unfortunately, the non-linear $u^4$ term takes away the strongest tool in our arsenal (Green's functions) so I don't have any clue how you'd proceed, or if an analytical solution is even possible (it likely isn't).
As one final note, I should say that I've made a few simplifying assumptions. Aside from the assumption of spherical symmetry, I've also implicitly assumed that there is only one radioactive species and that the density of rock is constant throughout the planet, among others. However, if you understand any of this and desire a more accurate approach, I should have hopefully given you enough to work with to extend it to a more accurate model.
Like I said before, this is a pretty complicated problem so don't be worried if you don't understand this approach-- I just wanted to leave a record for how this could be solved in an accurate manner. My apologies if it is also rather terse-- I don't have time right now to add more explanations. If I have more time later, I might edit some in.
