We can use grpdelay() to evaluate group delay of Discrete Time Filter in MATLAB.
However as far as I know, there is no such function for Continuous Time Filter in MATLAB.

I use following my custom function for this purpose without using ArcTan().

Do you have more cool method ?

function [Gdelay, fout] = my_Gdelay(H, fin)
x = real(H(:));
y = imag(H(:));
df = diff(fin(:));
dx = diff(x);
dy = diff(y);
x = x(1:end-1);
y = y(1:end-1);

Gdelay = -(dy.*x-y.*dx)./(x.^2+y.^2)./(2*pi*df);
fout = fin(1:end-1);