Taxicab numbers algorithm check With the programming language skills that are available to me at the time, I've written this program to find the "taxicab numbers" (e.g.  a number
expressible as the sum of two cubes in two different ways.)

While this code does work, it is definitely not scalable and it already takes about a minute to solve this for the below numbers. Does anyone have any tips/advice in order to make the complexity less?
Thanks in advance.

from math import pow

POWER = 3

for outcome in range(1000, 5000):
    maxtocheck = int( pow( outcome, 1/POWER ) )
    result = 0
    for i in range(0,maxtocheck+1):
        for j in range(0,maxtocheck+1):
            for k in range(0,maxtocheck+1):
                for l in range(0,maxtocheck+1):
                    if outcome == (pow(i, 3) + pow(j, 3)) == (pow(k, 3) + pow(l, 3)):
                        if i != k and i != l and j != k and j != l:
                            if result != outcome:
                                print("{} = {}**3 + {}**3 = {}**3 + {}**3".format(outcome,i,j,k,l))
                                result = outcome
As you said in the question you don't think this is very performant.

You only need two loops. If you go through every \$i\$ and \$j\$, then you have gone through every \$k\$ and \$l\$.
And so you don't need the \$k\$ and \$l\$ loops.

After this you have the equation:

$$\text{outcome} = i^3 + j^3$$

Since you know what \$\text{outcome}\$ is, it doesn't make sense for it to be what you're finding.
Instead if you re-arrange your code so that you can find \$j\$ for each \$i\$, then you can do everything in one loop.

$$j = (\text{outcome} - i^3)^{\frac{1}{3}}$$

And so to generate a list of \$i\$ and \$j\$ you can:

def generate_products(total):
    for i in range(int(total ** (1/3))+1):
        t = total - i ** 3
        # fix floating point error
        j = round(t ** (1/3))
        if t == j**3:
            yield i, j


After this you want to then combine these. If you're ok with \$i = k,j = l\$ then you can use itertools.product, otherwise you'd want itertools.permutations.
Which can produce all the numbers from 1 to 5000 in a couple of seconds.

def generate_products(total):
    for i in range(int(total ** (1/3))+1):
        t = total - i ** 3
        # fix floating point error
        j = round(t ** (1/3))
        if t == j**3:
            yield i, j

def generate_taxicab_numbers(outcome):
    for (i, j), (k, l) in itertools.permutations(generate_products(outcome), 2):
        print("{} = {}**3 + {}**3 = {}**3 + {}**3".format(outcome, i, j, k, l))

if __name__ == '__main__':
    for outcome in range(1000, 5000):
        generate_taxicab_numbers(outcome)


The slowest part of this is the print calls, and so if you move them out of the loops the code becomes even faster.

def generate_products(total):
    for i in range(int(total ** (1/3))+1):
        t = total - i ** 3
        # fix floating point error
        j = round(t ** (1/3))
        if t == j**3:
            yield i, j

def generate_taxicab_numbers(outcome):
    for (i, j), (k, l) in itertools.permutations(generate_products(outcome), 2):
        yield i, j, k, l

if __name__ == '__main__':
    to_print = []
    for outcome in range(1000, 5000):
        for nums in generate_taxicab_numbers(outcome):
            to_print.append("{} = {}**3 + {}**3 = {}**3 + {}**3".format(outcome, *nums))
    print('\n'.join(to_print))