A surprisingly hard-to-find answer. But I'm lucky that somebody else has asked almost the same question here.
I'm posting the answer that worked for me :
x = matrix(:,1);
y = matrix(:,2);
z = matrix(:,3);
xi=linspace(min(x),max(x),30)
yi=linspace(min(y),max(y),30)
[XI YI]=meshgrid(xi,yi);
ZI = griddata(x,y,z,XI,YI);
contourf(XI,YI,ZI)
which prints a nice color map.