Finding Ramanujan's taxicab numbers

I found the answer in Wolfram Mathworld

Taxicab[n_, max_] := {#[[1, 1]], First /@ Rest /@ #} & /@ 
Select[Split[
Sort[{Plus @@ #, #} & /@ Subsets[Range[Floor[max^(1/3)]]^3, {2}]],
First[#1] == First[#2] &], Length[#] == n &]

Taxicab[2, 10^5]

{{1729, {{1, 1728}, {729, 1000}}}, {4104, {{8, 4096}, {729, 3375}}}, {13832, {{8, 13824}, {5832, 8000}}}, {20683, {{1000, 19683}, {6859, 13824}}}, {32832, {{64, 32768}, {5832, 27000}}}, {39312, {{8, 39304}, {3375, 35937}}}, {40033, {{729, 39304}, {4096, 35937}}}, {46683, {{27, 46656}, {19683, 27000}}}, {64232, {{4913, 59319}, {17576, 46656}}}, {65728, {{1728, 64000}, {29791, 35937}}}}

First /@ %

${1729, 4104, 13832, 20683, 32832, 39312, 40033, 46683, 64232, 65728}$

ListPlot[First /@ Taxicab[2, 10^10], PlotStyle -> Red]

One way might be

n = 1729;
m = Floor[N[n^(1/3)]]
o[1] = Subsets[Range[m], {2}]
Cases[o[1], {x_, y_} /; x^3 + y^3 == n]

Another way

o[2] = Total[#^3] & /@ o[1]
Extract[o[1],Position[o[2], n]]

Here are the 2 taxi numbers found up to 10000

   o[3] = Last@Reap@Do[m = Floor[N[n^(1/3)]];
            o[1] = Subsets[Range[m], {2}];
            o[2] = Cases[o[1], {x_, y_} /; x^3 + y^3 == n];
            If[o[2] =!= {} && Length[o[2]] >= 2,
                 Sow[{n, o[2]}]
            ],
            {n, 1729, 10000, 1}
            ]

Grid[Flatten[o[3], 1], Frame -> All]

Verify