Rather than go through the arduous development of Taylor's theorem for functions of two variables, I'll say a few words and then present the theorem. In the one variable case, the nth term in the approximation is composed of the nth derivative of the function. For functions of two variables, there are n+1 different derivatives of nth order. For example, fxxxx, fxxxy, fxxyy, fxyyy, fyyyy are the five fourth order derivatives. There are actually more, but due to the equality of mixed partial derivatives, many of these are the same. Thus, our formula for Taylor's theorem must incorporate more than one derivative at each order.
The formula for a third order approximation to f(x,y) near (x0,y0) is
The factors of 2 and 3 appearing the second and third order mixed partial terms are due to the fact that there are two equal mixed partials derivatives of second order and a pair of three equal third order mixed partials.
When we start trying to maximize and minimize functions of two variables, we'll see that, near a maximum or minimum, all functions look more or less like a quadratic. Thus, the second order, or quadratic approximation will be useful: