David devised a scheme to find a polynomial rule for any number of input and output pairs of numbers you come up with, from your rule. [He has written this up and it is in the hands of the NCTM for publication].

Don spent over a month trying to understand what David did. After 55 years of teaching math, Don had never seen anything like this. He worked hard to try examples using David's method, then encouraged David to write this up. 

Don made up the rule y=x + x, to test David's method and figured out the approximations for each y-value to 8 digit accuracy, using Mathematica.

N[0 + 0, 8]

1

N[1 + 1, 8]

3.71828

N[2 + 2, 8]

9.38906

N[3 + 3, 8]

23.085537

Don's rule y=x + x, gives the following pairs of numbers: (x,y), (0,1), (1,1 + 1), (2,2 + 2), (3,3 + 3)

He used the command Expand on David's work (not shown here) and got the following function: 

1 + 2.93311x - 1.06036x2 + 0.845536 x3

 

Don used Mathematica below to show that when the x-value is substituted into the function (using .. /.x-> ), the y-value comes out as it should.

1 + 2.9331069780437167`x - 1.0603608348801465`x2 + 0.8455356852954754`x3/.x0

1

1 + 2.9331069780437167`x - 1.0603608348801465`x2 + 0.8455356852954754`x3/.x1

3.71828

1 + 2.9331069780437167`x - 1.0603608348801465`x2 + 0.8455356852954754`x3/.x2

9.38906

1 + 2.9331069780437167`x - 1.0603608348801465`x2 + 0.8455356852954754`x3/.x3

23.0855

Don graphed David's found function (a cubic polynomial) and Don's exponential function in Mathematica, from 0 to 8. The graphs show that from 0-3 they are exactly the same, then Don's exponential function jumps up faster when x>3. So David's method of finding a rule for the 4 pairs of numbers -or any number of pairs- really works, and the arithmetic in David's scheme was much, much, easier using Mathematica.

Plot[{1+2.933106978043716`x-1.0603608348801465`x2 + 0.845535685295475`x3,x+x},{x,0,8}]