For 3D data cubes of things like ocean temperature, you most often need to permute the first two dimensions when reading in from a netcdf. That would look like this for your data:
sst = permute(ncread(ncfile, 'thetao'),[2 1 3]);
To make a 2D map of this 3D data, you'll need to either pick a single time slice and plot it like this:
to plot the third map of SST data in the time series. Or perhaps you wish to plot the average SST map for all of your data. That would look like this:
which calculates the mean long the third dimension (time) and omits any missing data values if there are any.
A side note: As a matter of style, it's generally best if you can avoid using the variable names X and Y to mean longitude and latitude, because X and Y often refer to projected coordinates which are in meters, rather than degrees. In my own work, I typically use lowercase lon and lat for 1D arrays and uppercase Lon and Lat for 2D grids of longitdue and latitude. That would look like
[Lon,Lat] = meshgrid(lon,lat);