Increase performance for Within- and Between-Cluster Scatter (Three nested loops?)

I want to implement the following two equations in Matlab to calculate the Within-Cluster-Scatter and the Between-Cluster-Scatter, with:

d(p,p') = distance between two points, I use the squared eucleadian distance

Uc = Number of points in cluster c

UT = Number of all points from all clusters

Sc = Aggregate of all points from cluster c

This is the function I wrote for this, but it takes quite some time and I assume that it is not using all the capabilites that matlab offers. But I wasn't able to find a way to get rid of the nested loops. Can anybody offer suggestistions for improvment?

function [B, W, R]=withinAndBetweenScatter(data,Y)

    % data  = m x n matrix, with m observations and n variables
    %         m in the order of a 800 or so
    %         n something between 1 and 80
    % Y     = vector with the clusters the observations belonging to

    clases=unique(Y);           % GET VECTOR OF CLASSES
    N_clases=length(clases);    % HOW MANY CLASSES
    
    B = 0;  % Between Scatter
    W = 0;  % Within Scatter          

    for i=1:N_clases
        clasei = find(Y==clases(i));    % GET DATA FOR EACH CLASS
        xi = data(clasei,:);            % all data from class i (Sc)
        li = size(xi,1);                % Number of data points in class i (Uc)
        
        claseni = find(Y~=clases(i));    
        xo = data(claseni,:);       % all data that is outside class i
        lo  = size(xo,1);           % Number of data points outside class i (Uc-Ut)
        
        sumW = 0;   % Within Sum
        % Sum over cost for all points with each other within cluster
        for j = 1:li 
            for k = 1:li
                if (k ~= j)
                    % Squared euclidean distance as dissimilarity measure 
                    sumW = sumW + (norm(xi(j,:) - xi(k,:)))^2;
                end
            end
        end   
        W = W + (sumW / li^2);
   
        
        sumB = 0;   % Between Sum
        % Sum over cost for all points with point points outside this
        % cluster
        for j = 1:li
            for k = 1:lo
                sumB = sumB + (norm(xi(j,:) - xo(k,:)))^2;
            end
        end
        B = B + (sumB / (li*lo));     
        
    end
    
    W = W/2;
    B = B/2;
    
    % Ratio of between cluster scatter to within cluster scatter
    R = B/W;
end