The temperature of any object in space is the equilibrium of incomming energy and outgoing energy. In the case of light it's the equilibrium between the energy which is absorbed and the energy which is radiated.
One can assume a lambertian surface
https://en.wikipedia.org/wiki/Lambertian_reflectance to make the computation simpler. The surface of the moon for example is
not a lambertian, there are more complicated BDRF's to compute it more exactly. But this isn't necessary if the solution doesn't have to be very accurate.
How to compute the incomming energy from a distant light source? One way is to compute a exact solution by using a solid angle
https://en.wikipedia.org/wiki/Solid_angle which is spanned by the projected circle in the angular coordinate system from the origin of the point light and do the math based on that.
A simpler way is to assume that the light source is very far away, so that the directions of the light are assumed to be parallel on all points of the receiving body. We can then compute the ratio of the surface of a sphere (the radius here is the distance between the light source and the other body) and the area of the projected sphere of the receiver (which is a circle with the radius of the receiver).
areaRatio = areaOfEmitterSphere / areaOfProjectedReceiverSphere
areaOfEmitterSphere = 4π*pow(distance,2)
areaOfProjectedReceiverSphere = pi*pow(r, 2)
We can then compute the received energy as the product of the total power of the source and the ratio:
receivedPower = sourcePower*areaRatio
So now we got a way on how to compute the receivedPower which heats up the body (rate of heating depends on heat capacity of the body).
So now we are left of with the problem on how to compute the outgoing power.
The only way a body can emit energy in space is by black body radiation
https://en.wikipedia.org/wiki/Black-body_radiation .
We can then plug in the temperature of the body to get the power which is radiated.
...
Now we are left of by finding a way on how to solve this nonlinear relationship to get the temerature where the equilibrium is maintained. This could be solved by using root finding algorithms, or simply simulating for N steps with a given timeframe to simulate per step.