Extend the one-variable optimizer to a two-dimensional surface f(x,y)=x²+y². Return every point so convergence can be inspected axis by axis.

def gradient_descent_2d(x, y, learning_rate, steps):
    return [[step, x, y, loss], ...]

• record the starting point before updating
• update both coordinates using their partial derivatives
• record the joint loss after every update