function factor(_n) {
	n = _n
	a = J(14, 2, .)
	i = 0
	if (mod(n, 2)==0) {
		j = 0
		while (mod(n, 2)==0) {
			j++
			n = n/2
		}
		i++
		a[i,1] = 2
		a[i,2] = j
	}
	for (k=3; k*k<=n; k=k+2) {
		if (mod(n, k)==0) {
			j = 0
			while (mod(n, k)==0) {
				j++
				n = n/k
			}
			i++
			a[i,1] = k
			a[i,2] = j
		}
	}
	if (n>1) {
		i++
		a[i,1] = n
		a[i,2] = 1
	}
	return(a[1::i,.])
}

function sumdigits(_n) {
	n = _n
	for (s=0; n>0; n=floor(n/10)) s = s+mod(n,10)
	return(s)
}

function smith(n) {
	a = J(n, 1, .)
	i = 0
	for (j=2; j<=n; j++) {
		f = factor(j)
		m = rows(f)
		if (m>1 | f[1,2]>1) {
			s = 0
			for (k=1; k<=m; k++) s = s+sumdigits(f[k,1])*f[k,2]
			if (s==sumdigits(j)) a[++i] = j
		}
	}
	return(a[1::i])
}

a = smith(10000)
n = rows(a)
n
  376

a[1::10]'

         1     2     3     4     5     6     7     8     9    10
    +-------------------------------------------------------------+
  1 |    4    22    27    58    85    94   121   166   202   265  |
    +-------------------------------------------------------------+

a[n-9::n]'

          1      2      3      4      5      6      7      8      9     10
    +-----------------------------------------------------------------------+
  1 |  9843   9849   9861   9880   9895   9924   9942   9968   9975   9985  |
    +-----------------------------------------------------------------------+