I don't have the mapping toolbox, so I'm just relying on the documentation to write my answer.
For a start, you can replace the inner loop with:
AzimuthPoints(ii,ii+1:end)=azimuth('rh',lat(ii),long(ii),lat(ii+1:end),long(ii+1:end));
I believe the following may work to replace both loops (and the AzimuthPoints declaration)
[lat1, lat2] = meshgrid(lat);
[long1, long2] = meshgrid(long);
AzimuthPoints = azimuth('rh', lat1, long1, lat2, long2);
The downside of this method is that it's going to calculate the whole AzimuthPoints matrix, not just the upper triangle, so it may actually be slower.
You could do:
to replace the lower triangle and main diagonal with zero, on the assumption that calculating
azimuth('rh', 0, 0, 0, 0)
is faster.
Best Answer