How could I find the maximum closed contour line?

Since a saddle point is where the contour lines go from closed to open, the idea is to find the location of the saddle point and evaluate the value of the streamline function there. That value is the maximum closed-streamline value, given that you are encircling a local minimum.

The method here uses a numerical gradient and not an analytical one, so occasionally the location might end up on a corner of the plot or somewhere like that. In that case you would just have to narrow the search to the appropriate region. Decreasing the spacing in x and y gives a more accurate answer but a denser plot.

% create streamline function
xx = -10:.4:10;
yy = -10:.4:10;
[x y] = meshgrid(xx,yy);
x0 = -5.1;
y0 = -3.5;
x1 = 4;
y1 = 3.1;
S = @(x,y) -sqrt(1./(5+(x-x0).^2+2*(y-y0).^2) + 1./(5+2*(x-x1).^2+(y-y1).^2));
% find where gradient of S is approximately zero
[gx gy] = gradient(S(x,y));
norm2 = gx.^2+gy.^2;  % function to minimize
[a ind] = min(norm2,[],'all','linear')
[j k] = ind2sub([numel(xx) numel(yy)],ind)
x0 = xx(k);   % saddle point location
y0 = yy(j);
S(x0,y0)    % max streamline value
% plot
figure(1)
quiver(x,y,gy,-gx)
hold on
contour(x,y,S(x,y))
%colorbar
plot(x0,y0,'o')
grid on
hold off