/* need to read in X1, X2, X3 first */
scalar n1=rowsof(X1)
scalar n2=rowsof(X2)
scalar n3=rowsof(X3)

mat u1=J(n1,1,1)
mat u2=J(n2,1,1)
mat u3=J(n3,1,1)

/* mv stands for mean vector */
mat mv1 = 1/n1*u1'*X1
mat mv2 = 1/n2*u2'*X2
mat mv3 = 1/n3*u3'*X3

/* mm stands for mean matrix */
mat mm1=u1*mv1
mat mm2=u2*mv2
mat mm3=u3*mv3

/* ds stands for deviation score matrix */
mat ds1=X1-mm1
mat ds2=X2-mm2
mat ds3=X3-mm3

/* sscp stands for the deviation sums of squares and cross products */
mat sscp1=ds1'*ds1
mat sscp2=ds2'*ds2
mat sscp3=ds3'*ds3

mat X=X1\X2\X3
scalar n = rowsof(X)
mat u = J(n,1,1)
mat mv = u'*X
mat mm = (1/n)*u*mv
mat td = X-mm  /* deviation scores */

mat tsscp = td'*td           /* total sscp */
mat w = sscp1+sscp2+sscp3    /* within sscp */
mat b = tsscp-w              /* between sscp */

mat wb = syminv(w)*b         /* w inverse b */
mat list wb

mat l = cholesky(w)          /* Cholesky decomposition */
mat u = l'                   /* u equals l transpose   */
mat a = inv(l)*b*inv(u)      /* a is symmetric */
mat symeigen uv e = a        /* e will have the eigenvalues of W-1B */
mat v = inv(u)*uv            /* v will have the eigenvectors of W-1B */
mat list e
mat list v

e[1,3]
            e1          e2          e3
r1   .89198797   .00524207  -7.851e-17

v[3,3]
            e1          e2          e3
c1   .06410227   .01442654   -.0314958
c2  -.00186162   .06888879   .05943386
c3    .0537507  -.02620577   .01270798

/* trimmed eigenvalues and eigenvectors */
mat e = e[1, 1..2]
mat v = v[1..., 1..2]
mat list e
mat list v

e[1,2]
           e1         e2
r1  .89198797  .00524207

v[3,2]
            e1          e2
c1   .06410227   .01442654
c2  -.00186162   .06888879
c3    .0537507  -.02620577