Try the JSL code below (open a script windom and paste and hit run).
sample_Q = Covariance( [35 36 40, 38 37 40, 48 41 42, 50 46 54] );
pop_Q = Covariance( [35 36 40, 38 37 40, 48 41 42, 50 46 54] ) * 3 / 4;
Show( sample_Q, pop_Q );
You should get the results
sample_Q =
[ 54,25 31,33 37,33,
31,33 20,66 28,66,
37,33 28,66 45,33];
pop_Q =
[ 40,6875 23,5 28,
23,5 15,5 21,5,
28 21,5 34];
The latter equals your calculation. The difference between the two is a factor of 3/4 i.e. (N-1)/1. I think the equation you have used is actually not for the (unbiased) sample covariance but rather the population covariance which would be the right choice if the population mean is known. Try to use N-1 as denominator instead of N and you should end up with the same covariance estimate as JMPs formula.