Did you look at the convex hull? MATLAB has a function convhull(). The points farthest apart will be on the convex hull.

If you have more points that you want than are on the convex hull, you'll have to consider interior points. For example if you have 3 points on the vertex of the triangle and a cluster of points interior to them, the convex hull is the triangle. But if you want the 5 pairs that are farthest from each other, the convex hull can give you only 3 so to get the other 2 you need, you'd probably have to do an exhaustive search finding the distance of every point to every other point and then sort. Fortunately if you have only a handful of points like you do, then this should not take very long.

In fact I want to maximize the total interpoint distance (something like the total distance of all the edges of a Delaunay Triangulation (or maybe more simply just the total sum of the distances between all the points)). But I want to avoid the convex hull, because then points are only selected on the edges leaving a hole in the center (I want to have points chosen well spread in my plane (i.e. also in the center if better))

Went with a "simple" way - tried to maximize the sum of the distances between each point in the sub-set of "n" below (similar to way you described earlier):

 n = randn(14,2);
 k = 7;

 b = nchoosek(1:length(n),k);
 c = zeros(size(b,1),1);

 for i = 1 : size(b,1)

    subb = b(i,:);
    subn = n(subb,:);

    for j = 1 : size(b,2)-1
        for k = j+1 : size(b,2)
            c(i) = c(i) + sqrt((subn(j,1)-subn(k,1))^2 + ...
                               (subn(j,2)-subn(k,2))^2);
        end
    end

 end

 [d,e] = max(c);

 f = b(e,:);

 %figure
 plot(n(:,1),n(:,2),'bo',n(f,1),n(f,2),'rx');
 legend('Point Set',sprintf('"Best" %d Points',k));

Elige Grant, your solution is very interesting, but it has the tendency to chose points on the exterior of the set of points, even if they are close one to the others. (e.g. in the figure you show, among the points selected in the top left corner, I would have expected the point in the middle not to be seleceted, but instead to have a point selected somewhere in the middle). So maybe just direct total sum of all the distances is not enough, a selection of distances like using a Delaunay triangulation might be better...

Based on Elige Grant's solution, I tried another solution using the Delaunay triangulation. It seems a bit better. (changing max_dist_crit to 1 and delny_tri_crit and max_mindist_crit to 0 allows to see Elige Grant's solution: points are selected from the outside to the inside, not the behaviour that I desired). But it was still not perfect. I like the solution below: maximizing the sum of minimum distance to the closest point for each point.

MATLAB code
% IF MAX LENGTH OF SEGMENTS OF DELAUNAY TRIANGULATION CRITERION:
delny_tri_crit=0;
% IF MAX OF ALL DISTANCES BETWEEN ALL POINTS WITHIN THE CHROMATOGRAM:
max_dist_crit=0;
% IF MAX LENGTH OF SEGMENTS OF DELAUNAY TRIANGULATION CRITERION, with
% double count of distances on the edge of the convex hull. Doesn't really
% work, even with other criteria than *2...
delny_tri_crit2=0;
% IF MAX OF the minimum DISTANCE to a neighboring peak for each alignment
% point:
max_mindist_crit=1;
Peaks_Ref=[ 98 1547; % (x,y) position
641 1108 ; 124 476 ; 508 507 ; 619 512 ; 746 531 ; 342 507 ; 439 1018 ; 195 1099 ; 550 843; 721 1651; 384 1547; 305 2364 ; 175 1649 ; ];
%TRI = delaunay(X,Y)
n = Peaks_Ref;
isfirst=true;
for k=3:size(Peaks_Ref,1)-1
%k = 11;
b = nchoosek(1:length(n),k);
c = zeros(size(b,1),1);
if max_mindist_crit
  c=c+1e99;
end
for i = 1 : size(b,1)
subb = b(i,:);
subn = n(subb,:);
if or(delny_tri_crit,delny_tri_crit2)
TRI = delaunay(subn(:,1),subn(:,2));
for ka=1:size(TRI,1)
for k2=1:3
c(i)=c(i)+sqrt(((subn(TRI(ka,k2),1)-subn(TRI(ka,mod(k2,3)+1),1))^2)+((subn(TRI(ka,k2),2)-subn(TRI(ka,mod(k2,3)+1),2))^2));
end
end
CVHL=convhull(subn(:,1),subn(:,2));
for kx=1:length(CVHL)-1
  if delny_tri_crit
c(i)=c(i)+sqrt(((subn(CVHL(kx),1)-subn(CVHL(kx+1),1))^2)+((subn(CVHL(kx),2)-subn(CVHL(kx+1),2))^2));
  elseif delny_tri_crit2
c(i)=c(i)+3*sqrt(((subn(CVHL(kx),1)-subn(CVHL(kx+1),1))^2)+((subn(CVHL(kx),2)-subn(CVHL(kx+1),2))^2));
  end
end
elseif max_dist_crit
for j = 1 : size(b,2)-1
for ka = j+1 : size(b,2)
    c(i) = c(i) + sqrt((subn(j,1)-subn(ka,1))^2 + ...
      (subn(j,2)-subn(ka,2))^2);
end
end
elseif max_mindist_crit
for j = 1 : size(b,2)-1
for ka = j+1 : size(b,2)
    c(i) = min(c(i),sqrt((subn(j,1)-subn(ka,1))^2 + ...
      (subn(j,2)-subn(ka,2))^2));
end
end 
end
end
[d,e] = max(c);
f = b(e,:);
figure
plot(n(:,1),n(:,2),'bo',n(f,1),n(f,2),'rx');
legend('Point Set',sprintf('"Best" %d Points',k))
%on ajoute l'image à la fin de notre gif
F = getframe(gcf); %on récupère l'image
[RGB,badmap] = frame2im(F); %on la convertit en image de type 'true-color'
[IND,map] = rgb2ind(RGB, 255); %on convertit en couleur indéxées. 255 est le nombre de couleur.
if isfirst
  if delny_tri_crit
imwrite(IND,map,'FarthestAlPtsDelaunay.gif','gif'); %c'est la première image: on crée le fichier
  elseif max_dist_crit
   imwrite(IND,map,'FarthestAlPtsMaxDist.gif','gif'); %c'est la première image: on crée le fichier   
  elseif delny_tri_crit2
     imwrite(IND,map,'FarthestAlPtsDelaunay2.gif','gif'); %c'est la première image: on crée le fichier
 elseif max_mindist_crit
     imwrite(IND,map,'FarthestAlPtsMaxMin.gif','gif'); %c'est la première image: on crée le fichier
  end
isfirst=false;
else
%%%%%imwrite(IND,map,'example.gif','gif','WriteMode','append'); %ce sont les suivantes: on les ajoute
if delny_tri_crit
imwrite(IND,map,'FarthestAlPtsDelaunay.gif','gif','WriteMode','append','DelayTime',1); 
elseif max_dist_crit
imwrite(IND,map,'FarthestAlPtsMaxDist.gif','gif','WriteMode','append','DelayTime',1);
elseif delny_tri_crit2
imwrite(IND,map,'FarthestAlPtsDelaunay2.gif','gif','WriteMode','append','DelayTime',1); 
elseif max_mindist_crit
imwrite(IND,map,'FarthestAlPtsMaxMin.gif','gif','WriteMode','append','DelayTime',1); 
end
end
close
end

Added a few more tests (including your latest) and added your data (instead of just random dataset), which may be of interest. Relevant code:

 n = [ 98 1547; 641 1108 ; 124 476 ; 508 507 ; 619 512 ; 746 531 ; 342 507 ; 439 1018 ; 195 1099 ; 550 843; 721 1651; 384 1547; 305 2364 ; 175 1649]
 k = 7;

 b = nchoosek(1:length(n),k);
 c1 = zeros(size(b,1),1);
 c2 = zeros(size(b,1),1);
 c3 = zeros(size(b,1),1);
 c4 = zeros(size(b,1),1);
 c5 = zeros(size(b,1),1);

 for i = 1 : size(b,1)

    subb = b(i,:);
    subn = n(subb,:);
    subd = delaunay(subn(:,1),subn(:,2));

    % Method 1 (Sum of distances between each pair in points sub-set)
    for j = 1 : size(b,2)-1
        for k = j+1 : size(b,2)
            c1(i) = c1(i) + sqrt((subn(j,1)-subn(k,1))^2 + ...
                                 (subn(j,2)-subn(k,2))^2);
        end
    end

    % Method 2 (Sum of Delaunay triangle perimeter, squared)
    for j = 1 : size(subd,1)
        c2(i) = c2(i) + ... 
                ( sqrt((subn(subd(j,1),1)-subn(subd(j,2),1))^2 ...
                     + (subn(subd(j,1),2)-subn(subd(j,2),2))^2) ...
                + sqrt((subn(subd(j,3),1)-subn(subd(j,2),1))^2 ...
                     + (subn(subd(j,3),2)-subn(subd(j,2),2))^2) ...
                + sqrt((subn(subd(j,1),1)-subn(subd(j,3),1))^2 ...
                     + (subn(subd(j,1),2)-subn(subd(j,3),2))^2) )^2;
    end

    % Method 3 (Sum of Delaunay triangle area)
    for j = 1 : size(subd,1) 
        c3(i) = c3(i) + polyarea(subn(subd(j,:),1),subn(subd(j,:),2));

    end

    % Method 5 (Delaunay + Convhull)
    for ka=1:size(subd,1)
        for k2=1:3
            c5(i)= c5(i) + ...
            sqrt(((subn(subd(ka,k2),1)-subn(subd(ka,mod(k2,3)+1),1))^2)+...
                 ((subn(subd(ka,k2),2)-subn(subd(ka,mod(k2,3)+1),2))^2));
        end
    end
    CVHL=convhull(subn(:,1),subn(:,2));
    for kx=1:length(CVHL)-1
        c5(i)=c5(i) + ...
            sqrt(((subn(CVHL(kx),1)-subn(CVHL(kx+1),1))^2)+...
                 ((subn(CVHL(kx),2)-subn(CVHL(kx+1),2))^2));
    end

 end

 % Method 4 (Method 2 and Method 3)
 c4 = c2/max(c2) + c3/max(c3);

 [~,e1] = max(c1);
 [~,e2] = max(c2);
 [~,e3] = max(c3);
 [~,e4] = max(c4);
 [~,e5] = max(c5);

 f1 = b(e1,:);
 f2 = b(e2,:);
 f3 = b(e3,:);
 f4 = b(e4,:);
 f5 = b(e5,:);

 figure;
subplot(2,3,1);
    plot(n( :,1),n(: ,2),'bo','MarkerSize',15);
    title('Starting Point Set'); axis equal; 
subplot(2,3,2);
    plot(n( :,1),n(: ,2),'bo','MarkerSize',15); hold on;
    plot(n(f1,1),n(f1,2),'rx','MarkerSize',15); hold off;
    title(sprintf('"Best" %d Points (M-1)',k)); axis equal; 
subplot(2,3,3);
    plot(n( :,1),n(: ,2),'bo','MarkerSize',15); hold on;
    plot(n(f2,1),n(f2,2),'rx','MarkerSize',15); hold off;
    title(sprintf('"Best" %d Points (M-2)',k)); axis equal; 
subplot(2,3,4);
    plot(n( :,1),n(: ,2),'bo','MarkerSize',15); hold on;
    plot(n(f3,1),n(f3,2),'rx','MarkerSize',15); hold off;
    title(sprintf('"Best" %d Points (M-3)',k)); axis equal; 
subplot(2,3,5);
    plot(n( :,1),n(: ,2),'bo','MarkerSize',15); hold on;
    plot(n(f4,1),n(f4,2),'rx','MarkerSize',15); hold off;
    title(sprintf('"Best" %d Points (M-4)',k)); axis equal; 
subplot(2,3,6);
    plot(n( :,1),n(: ,2),'bo','MarkerSize',15); hold on;
    plot(n(f5,1),n(f5,2),'rx','MarkerSize',15); hold off;
    title(sprintf('"Best" %d Points (M-5)',k)); axis equal;