Matlab: Error while trying to extract data from a netCDF for a specific range

Question

I want to read the variable “tref” for desired latitude and longitude (which I index), 24th M, 4th L, and all S from a NetCDF file:

filename = 'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NCEP-CFSv2/.HINDCAST/.MONTHLY/tref/1+index/.tref/tref/tref/dods';
ncdisp(filename) % so the structure of tref is (X*Y*M*L*S which means lon*lat*ensemble*leadtime*date) 
lon = ncread(filename, 'X');
lat = ncread(filename, 'Y');
date = ncread(filename, 'S');
leadtimes = ncread(filename, 'L');
ensembels = ncread(filename, 'M');
ind_lon = find(lon>=44 & lon<=64);
ind_lat = find(lat>=24 & lat<=40);
tref = ncread(filename, 'tref', [ind_lon ind_lat 24 4 inf]);

I used the following command, but it was not successful. The error is:

Error using horzcat Dimensions of arrays being concatenated are not consistent.

Here are the results of the ncdisp()

Source:
           http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NCEP-CFSv2/.HINDCAST/.MONTHLY/tref/1+index/.tref/tref/tref/dods
Format:
           classic
Global Attributes:
           Conventions = 'IRIDL'
Dimensions:
           L = 10
           M = 24
           S = 348
           X = 360
           Y = 181
Variables:
    M   
           Size:       24x1
           Dimensions: M
           Datatype:   single
           Attributes:
                       standard_name = 'realization'
                       long_name     = 'Ensemble Member'
                       pointwidth    = 1
                       gridtype      = 0
                       units         = 'unitless'
    X   
           Size:       360x1
           Dimensions: X
           Datatype:   single
           Attributes:
                       standard_name = 'longitude'
                       pointwidth    = 1
                       gridtype      = 1
                       units         = 'degree_east'
    L   
           Size:       10x1
           Dimensions: L
           Datatype:   single
           Attributes:
                       standard_name = 'forecast_period'
                       long_name     = 'Lead'
                       pointwidth    = 1
                       gridtype      = 0
                       units         = 'months'
    S   
           Size:       348x1
           Dimensions: S
           Datatype:   single
           Attributes:
                       standard_name = 'forecast_reference_time'
                       long_name     = 'Forecast Start Time'
                       pointwidth    = 0
                       calendar      = '360'
                       gridtype      = 0
                       units         = 'months since 1960-01-01'
    Y   
           Size:       181x1
           Dimensions: Y
           Datatype:   single
           Attributes:
                       standard_name = 'latitude'
                       pointwidth    = 1
                       gridtype      = 0
                       units         = 'degree_north'
    tref
           Size:       360x181x24x10x348
           Dimensions: X,Y,M,L,S
           Datatype:   single
           Attributes:
                       pointwidth    = 0
                       standard_name = 'air_temperature'
                       PTVersion     = 2
                       long_name     = 'Reference Temperature'
                       process       = 'Spectral Statistical Interpolation (SSI) analysis from "Final" run.'
                       PDS_TimeRange = 3
                       gribNumBits   = 21
                       center        = 'US Weather Service - National Met. Center'
                       grib_name     = 'TMP'
                       gribcenter    = 7
                       scale_min     = 206.4713
                       gribparam     = 11
                       scale_max     = 315.4577
                       GRIBgridcode  = 3
                       gribleveltype = 105
                       gribfield     = 1
                       units         = 'Kelvin_scale'

Can you help me with this issue, please?