If I have the statistical uncertainties of the ensemble average magnetisation and the average energy from a monte carlo simulation of an Ising Model, how do I find the errors on the specific heat capacity and the magnetic susceptibility?
Assuming that you have magnetization as a function of the external field at constant temperature and energy as a function of temperature at constant external field, and that your Monte Carlo is a standard Monte Carlo where temperature and external field are fixed parameters, for both cases the problem is the estimate of error for numerical differentiation of noisy data. If you look at the literature, with the above keywords, you'll realize that the problem is ill-posed without specifying which method is used for the numerical differentiation.
Finite difference methods are quite dangerous when used with exact data and become almost useless with noisy data. Safer procedures involve some smoothing of the data like fitting them with a smoothing spline or by Fourier noise filtering. Then the derivative can be taken on the smooth approximation and the error analysis of the derivative becomes equivalent to error propagation from the errors on the smooth representation.
Unfortunately, it is difficult to provide explicit formulae without additional information on the planned smoothing procedure which is somewhat depending on the quality of data and on the required level of accuracy. However, some packages for smoothing spline fitting could provide error estimates (although I do not have immediately available references).