Solution (source code)

= Solution

The <partial derivatives> are $f_x=e^{xy}(1+xy)$ and $f_y=x^2e^{xy}$. At $(1,1)$ they have values $2e$ and $e$. The first-order <Multivariate Taylor polynomial> therefore gives the change
$$
\boxed{\delta f=2e\,\delta x+e\,\delta y+O(\delta x^2+\delta y^2).}
$$
Thus the requested linear approximation is $e(2\delta x+\delta y)$; it is the change from $f(1,1)=e$, not the new <function> value itself.