I am trying to compute a 1000 x 60 matrix or list of lists (and ideally this should go up to 1000 x 500 or 1000 x 1000).

Each element is the result of a FindRoot operation, so I make my list by doing

Table[Flatten[h /. FindRoot[h == F[h, b, g], b, 1, 1000, g, 1, 60)

but 16GB of RAM are filled up. I think I should be able to hold list of lists much bigger than that, so probably using Table with FindRoot is causing Mathematica to store a lot of undeeded stuff in memory.

Here is the code:

ι[m_, n_] := Binomial[n, n*(1 - m)/2]*2^(-n);
f[m_, h_, b_, g_, n_] := (h*m + g/2*m^2) + 
 1/(n*b)*Log[ι[m, n]];
μ[m_, h_, b_, g_, n_] := 
 Exp[b*n*f[m, h, b, g, n] + b*n*(-h + g/2)]/
 Sum[Exp[b*n*f[x, h, b, g, n] + b*n*(-h + g/2)], x, -1 + 2/n, 
 1 - 2/n, 2/n];
moment[h_, x_, b_, g_, n_] := Sum[m^x*μ[m, h, b, g, n], m, -1 + 2/n, 1 - 2/n, 2/n];
var[h_, b_, g_, n_] := moment[h, 2, b, g, n] - moment[h, 1, b, g, n]^2;
cov[h_, b_, g_, n_] := moment[h, 3, b, g, n] - moment[h, 1, b, g,n]*moment[h, 2, b, g, n];
F[h_,b_,g_,n_]:= -d*b*(cov[h, b, gg, n] + 
 2 var[h, b, gg, n]);
n = 100;
d = 0.9;

glist = Table[g, g, 0.4, 1, 0.01];
blist = Table[b, b, 1.1, 10.1, 0.01];

heatdata = Table[
 Flatten[h /. 
 FindRoot[
 h == F[h,b,g,n], h, -0.01]][[1]]
 , b, blist, g, glist];

Your function F is implemented really, really inefficiently. By quite simple means and in the proposed situation, it can be sped up by a factor of 20000. The key is to start with calculations in machine precision as early as possible and to store frequently used data in packed arrays.

n = 100;
mlist = Range[-1. + 2/n, 1. - 2/n, 2./n];
m2list = mlist^2;
m3list = mlist^3;
logiotalist = Log[Binomial[n, n*(1 - mlist)/2]*2^(-n)];

d = 0.9;
glist = Range[0.4, 1, 0.01];
blist = Range[1.1, 10.1, 0.01];

ClearAll[F];
F[h_?NumericQ, b_, g_] := 
 Module[var, cov, explist, μlist, mom1, mom2, mom3,
 explist = Exp[(b n h) mlist + (b n g/2) m2list + logiotalist + b n (-h + g/2)];
 μlist = explist/Total[explist];
 mom1 = μlist.mlist;
 mom2 = μlist.m2list;
 mom3 = μlist.m3list;
 var = Subtract[mom2, mom1 mom1];
 cov = Subtract[mom3, mom1 mom2];
 (-d b) (cov + 2. var)
 ];

Just a quick test for precision and speed:

t1, r1 = F[0.1, blist[[1]], glist[[1]], n] // RepeatedTiming;
t2, r2 = Fnew[0.1, blist[[1]], glist[[1]]] // RepeatedTiming;
Abs[r1 - r2]/r1
t1/t2

-1.32375*10^-14

2.1*10^4

Now the parallelized solve loop requires about 10 seconds on my Quad Core Haswell CPU:

ParallelEvaluate[Off[General::munfl]];
heatdata = Developer`ToPackedArray[
 ParallelTable[
 Block[h0, h,
 h0 = -0.01;
 Developer`ToPackedArray[
 Table[
 h0 = h /. FindRoot[h == Fnew[h, b, g], h, h0],
 b, blist]
 ]
 ],
 g, glist]
 ]; // AbsoluteTiming // First

10.072

Memory considerations

You also see: Limited amount of RAM is not an issue here. That must have been caused by excessive memoziation. For the timing, it is crucial how information is stored and retrieved.
Storing computed values in a packed array for retrieving them later is significantly more efficient than memoization. Memoization into DownValues uses a complex data structure such as a hash table at its backend, and this data structure has certain overhead. In contrast, a packed array represents basically a connected block of physical memory, accompanied by some bytes of meta information (array dimensions and maybe some row pointers). Moreover, computation with data stored in packed arrays can take advantage of vectorization, which is most crucially employed in the following line:

explist = Exp[(b n h) mlist + (b n g/2) m2list + logiotalist + b n (-h + g/2)];

Remark on precision

Finally, I have to note that there is numerical underflow occurring in the course of the computation. This is probably caused by calling Exp with negative numbers of oversized absolute value. I decided to turn off the warning message, but this may lead to a significant loss of precision. So use with care. If one wants to do it correctly, one should investigate this further and apply, e.g. Clip or Threshold.