I was running into problems when using my fitted terms in another program: the same coefficients did not give the same surface.
I found the cause: the cosine and sine terms are inverted in zernfun. M>0 should be a cosine, while M<0 should be a sine.

I've corrected this, starting at line 193

if any(idx_pos)
z(:,idx_pos) = y(:,idx_pos).*cos(theta*m(idx_pos)');
end
if any(idx_neg)
z(:,idx_neg) = y(:,idx_neg).*sin(theta*-m(idx_neg)');
end

% note the required sign change "-m" in the sine term