replicate Robust Standard Errors with formula

Hello, I am trying to manually calculate robust SEs for a heteroskadistic data set.  
All of the formulas I have attempted generate either a 500X1 or 500X500 matrix. 
How can I compare my results with the 2X1 matrix generated by WhiteSE command? 

Please see below my code.  All suggestions are welcomed.

Thanks in advance, 
Robert 

    clear
    n=500;
    X=[1:500];
    y=2*X+normrnd(0,sqrt(X.^1.8));
    mdl=fitlm(X,y,'linear'); %add regression line
    b=X\y; % OLS estimate
    e=y-X*b; % OLS residuals
    SSRb=e'*e; % SSR function at OLS estimate
    %plot(X,y,'o')
    plot(mdl)
    plotResiduals(mdl,'fitted')
    
    %%Unadjusted statistics 
    h_cov = mdl.CoefficientCovariance;
    h_SE = mdl.SSE;
    h_coeff = mdl.Coefficients;
    
    %%Robust Standard Errors
    [WhiteCov,WhiteSe,Whitecoeff] = hc(mdl); % White's estimates
    
    %HCE= 1/n*(SSRb/e).*(SSRb/e.*(e.^2)).*(SSRb/e); 
    %HCE= 1/n*(sum((SSRb/e)).*(sum((SSRb/e.*(e.^2))).*(sum(SSRb/e)); 
    %HCE= 1/n*(sum((SSRb/e)).*(sum((SSRb/e.*(e.^2))).*(sum(SSRb/e))));
    %HCE = (SSRb/e).*SSRb'.*(diag(SSRb)).*(SSRb/e);
    %HCE = (SSRb/e).*SSRb'.*(diag(e.^2)).*(SSRb/e);
    HCE = (SSRb/-1).*(SSRb'.*(diag(e.^2))).*(SSRb/-1);