I am trying to use Taylor series to approximate the value of $f(x,y) = \frac{1}{x+y}$ for some code I am writing. Using WolframAlpha to quickly compute the Maclaurin series gives me: $$ \frac{1}{y} - \frac{x}{y^2} + \frac{x^2}{y^3} - \frac{x^3}{y^4} + ... $$
However, expanding $f(x,y)$ by hand using the following: $$ f(x,y)=f(x_0,y_0)+[f'_x(x_0,y_0)(x-x_0)+f'_y(x_0,y_0)(y-y_0)]+\frac{1}{2!}[f''_{xx}(x_0,y_0)(x-x_0)+2f''_{yx}(x_0,y_0)(x-x_0)(y-y_0)+f''_{yy}(x_0,y_0)(y-y_0)^2]+... $$
I obtained this result (substituting $x_0 = 0$): $$ \frac{1}{y_0} - \left[\frac{x}{y_0^2} + \frac{y-y_0}{y_0^2}\right] + \frac{1}{2!}\left[\frac{2x^2}{y_0^3} + 2\frac{2x(y-y_0)}{y_0^3} + \frac{2(y-y_0)^2}{y_0^3} \right] + ... $$
It appears WolframAlpha made the assumption $y=y_0$. I am wondering why/when is this a valid assumption to make and how does it affect the approximation of $f(x,y)$, relative to using the Maclaurin series computed by hand.
A small remark : it is not a Maclaurin series but a Taylor series and this may explain part of your problems.
If you have a look here, you will notice that your formula is correct. Around a point $(a,b)$, as given in the linked page, we have $$f(x,y)=f(a,b) +(x-a) f_x(a,b) +(y-b) f_y(a,b) + $$ $$\frac{1}{2!}\Big( (x-a)^2 f_{xx}(a,b) + 2(x-a)(y-b) f_{xy}(a,b) +(y-b)^2 f_{yy}(a,b) \Big)+\cdots$$ So, making $a=0$, you should get $$f(x,y)=f(0,b) +x f_x(0,b) +(y-b) f_y(0,b) + $$ $$\frac{1}{2!}\Big( x^2 f_{xx}(0,b) + 2x(y-b) f_{xy}(0,b) +(y-b)^2 f_{yy}(0,b) \Big)+\cdots$$ For your function $$f(0,b)=\frac{1}{b} \qquad f_x(0,b)=f_y(0,b)=-\frac{1}{b^2}\qquad f_{xx}(0,b)=f_{xy}(0,b)=f_{yy}(0,b)=\frac{2}{b^3}$$ which make $$\frac{1}{x+y}=\frac{1}{b}-\frac{1}{b^2}\left(x+(y-b)) \right)+\frac{1}{b^3}\left(x^2+2x(y-b)+(y-b)^2\right)+\cdots$$ which is exacly what you properly wrote.
May I confess that I am reluctant to use CAS when facing expansions in several variables ?