Plotting the time evolution of systems of differential equations of motion

Notice first that we can write each differential (dy/dx)_n as follows:
suppose y, y1, y2, y3,... are the last n values the function takes on over a certain N x values separated by a delta value.
A first derivative can be approximated like:

A second like:

And so on using divided differences of the function over a certain range.
So just picking a random relationship in x and y versus t might be:

Plugging in the approximations from above gives:

we can solve both of these for x and y, respectively, then plug in our delta for t, which I will make 1/10.:

Now we see we need to know the last two values of x (x1 and x2) and y (y1 and y2) to figure out the new x and y value.
Let's suppose the function started x1 = .5, x2 = .4, y1 = .75, y2 = 1
x = 0.5721630376
y = .7621630376
Then make these two values our x1 and y1, make the old x1 and y1: x2 and y2, and iterate. A plot of this looks like: