But nowhere either there or elsewhere online can I find how to calculate the profile drag.
Seems like you were searching for the incorrect terms. Try searching "boundary layer integral momentum equation", "von Karman boundary layer integral", "momentum deficit", "boundary layer Kinetic energy dissipation integral" or something similar. You should get plenty of resources/lecture notes on the subject.
So does anybody know how to calculate the form drag for a given aerofoil (eg NACA4412)?
In essence, this is done via integrating the "Integral Momentum Equation" or its non-dimensionalized version "von Karman Equation" over a coordinate system defined over the airfoil contour.
I can get the skin friction drag from boundary layer calculations but have no idea where to go from there.
If you study the "Momentum integral equation" stated above you can see there are two contributors to the momentum deficit or its rate of change. one is skin friction and the other is the profile drag which is equal to - mass defect * rate of change of edge velocity.
Put it in another way, the skin friction term usually is dominated by the front portion of an airfoil where due/ds mostly favorable and the profile drag is dominated by the pressure recovery region aft of the airfoil because of the negative velocity gradient and higher mass defect.
But the story does not end there, because you will have to add the viscous dissipation happening in the wake up to the far-field as well. This is usually done using the "Squire-Young" equation by applying the previously obtained momentum deficit at the trailing edge.
All above being said, you can easily use XFOIL for such calculations. But understanding what's going on under the hood in such an application, might provide one with great insight and appreciation for what they have got for free.