nasa = wmsfind('NASA Earth Observations','SearchField','any');
modis = refine(nasa,'land*day*month');
modis = modis(1);

server = WebMapServer(modis.ServerURL);

updatedLayer = updateLayers(server,modis);

mapRequest = WMSMapRequest(updatedLayer,server);
A = getMap(server,mapRequest.RequestURL);
R = mapRequest.RasterReference;

figure
ax = worldmap('world');
geoshow(A,R)
setm(ax,'MLabelParallel',-90,'MLabelLocation',90)
title({'MODIS Global Mosaic',modis.LayerTitle})

usgsLayers = wmsfind('usgs*imagery','SearchField','serverurl');

updatedLayers = [];
serverURLs = usgsLayers.servers;

for k=1:numel(serverURLs)
   serverLayers = refine(usgsLayers, serverURLs{k}, ...
          'SearchField','serverurl','MatchType','exact');
   serverURL = serverLayers(1).ServerURL;
   server = WebMapServer(serverURL);
   try
      layers = updateLayers(server,serverLayers);
      % Grow using concatenation because layers can have any
      % length ranging from 0 to numel(serverLayers).
      updatedLayers = [updatedLayers; layers];
   catch e
       fprintf('Server %s is not responding.\n', ...
              serverURL);
       fprintf('Error message is %s\n', e.message)
    end
end