simple_distance.gnuplot

set terminal pngcairo font "arial,8" size 640,350
set output 'simple_distance.png'
set multiplot layout 1, 2

r(d) = d*pi/180
dist_hav(x,y,x0,y0) = 2*asin(sqrt(sin(r(y-y0)/2)**2 + cos(r(y0))*cos(r(y))*sin(r(x-x0)/2)**2))*6371
dist_sim(x,y,x0,y0) = sqrt(r(y-y0)**2 + (cos(r(y0))*r(x-x0))**2)*6371
error(x,y,x0,y0) = 100*abs(dist_hav(x,y,x0,y0)-dist_sim(x,y,x0,y0))/dist_hav(x,y,x0,y0)

set view map 
set palette cubehelix negative
set contour
set isosamples 200, 200 
set cntrparam level discrete 5

set title "Relative error at 89°N"
set xrange [-.05/cos(r(89)): .05/cos(r(89))]
set yrange [89-.05:89+.05]
splot error(x,y,0,89) notitle with pm3d, \
    dist_hav(x,y,0,89) nosurface notitle

set title "Relative error at 65°N"
set xrange [-.05/cos(r(65)):.05/cos(r(65))]
set yrange [64.95:65.05]
splot error(x,y,0,65) with pm3d notitle, \
    dist_hav(x,y,0,65) nosurface notitle

unset multiplot

This is awesome! You can become even more accurate if you use the formula for equirectangular distance calculation from https://www.movable-type.co.uk/scripts/latlong.html :

earth_radius = 6371
dist_eqi_x(lng0, lat0, lng1, lat1) = r(lng1 - lng0) * cos(r(lat1 + lat0)/2)
dist_eqi_y(lat0, lat1) = r(lat1 - lat0)
dist_eqi(lng0, lat0, lng1, lat1) = earth_radius * sqrt(dist_eqi_x(lng0, lat0, lng1, lat1)**2 + dist_eqi_y(lat0, lat1)**2)