wlog, let a ≤ b ≤ c

For the subset of solutions where a=b=c, the sod(a) must be the square root of a and thus a is a perfect square.   A quick search for n up to 1 million shows that only 1 and 9 have this property.  sod(81) = 9.   For a d-digit number, its square will have about 2d digits with a maximum sod of 18d, but the d-digit number is on the order of 10^d or 10^(d-1).  Setting 18d = 10^d has a solution of d = approx 1.4, and 10^1.4 = approx 25.  So no n>25 can have the sod of its square greater than or equal to n.  In fact, 17 is the largest n where this happens:  sod(289) = 19 > 17.

So we have two solutions where a=b=c:

(1,1,1) and (81,81,81)

I suspect a similar argument can be made for an upper limit on the size of a,b,c in the general case where they are not all equal.

Suppose c is the largest and that it is a d-digit number.

c = sod(a)*sod(b)

Both a and b have a maximum of d digits.

The largest sod(a) can be is 9d.

The upper limit of c is 81*d^2

Can c be a 4-digit number?  (4*9)^2 = 1296 which is 4 digits.

If c could be a 5-digit number, then the largest a and b could be would be 99999 each, each with sod of 45; 45^2 is 2025 which is only 4 digits.

This sets an upper limit on d of 4 digits and on c of 81*4*4 = 1296

A computer search up to 1296 showed one more.

---------------

def sod(n):

    """ Input an integer.  Returns the Sum of the Digits  """

    aList = list(str(n))

    ans = 0

    for c in aList:

        ans = ans + int(c)

    return ans

big = 1296

for a in range(1,big+1):

    for b in range(a,big+1):

        for c in range(b,big+1):

            sa = sod(a)

            sb = sod(b)

            sc = sod(c)

            if a == sb*sc and b == sa*sc and c == sa*sb:

                print('({},{},{})'.format(a,b,c))