Formulation of -div(k*grad u)=f in Omega for the Finite Element Method

% Example1a
%
%   This script demonstrates the use of the Fem code to
%   solve the BVP
%
%          -div(k*grad u)=f in Omega,
%                       u=g on Gamma_1,
%                 k*du/dn=h on Gamma_2.
%
%   Here Omega is the unit square, Gamma_1 consists of
%   the top and left edges of the square and Gamma_2 is
%   the remainder of the boundary (the bottom and right
%   edges).  The coefficient is k(x,y)=1+x^2*y, and
%   f, g, h are chosen so that the exact solution is
%   u(x,y)=exp(2*x)*(x^2+y^2).

%   This routine is part of the MATLAB Fem code that
%   accompanies "Understanding and Implementing the Finite
%   Element Method" by Mark S. Gockenbach (copyright SIAM 2006).

% Define the problem functions (note that u and its partial derivatives
% are used here to compute the boundary values; also, these functions
% are used to check the accuracy of the computed solution):

f=vectorize(inline('-2*x*y*(2*exp(2*x)*(x^2+y^2)+2*exp(2*x)*x)-(1+x^2*y)*(4*exp(2*x)*(x^2+y^2)+8*exp(2*x)*x+2*exp(2*x))-2*x^2*exp(2*x)*y-2*(1+x^2*y)*exp(2*x)','x','y'));

This is the first example of some MatLab code of Mark Gockenbach’s. He declares that -div(k*grad u)=f. In the first line of MatLab code he declares

f=vectorize(inline('-2*x*y*(2*exp(2*x)*(x^2+y^2)+2*exp(2*x)*x)-(1+x^2*y)*(4*exp(2*x)*(x^2+y^2)+8*exp(2*x)*x+2*exp(2*x))-2*x^2*exp(2*x)*y-2*(1+x^2*y)*exp(2*x)','x','y'));

So there must be some equivalency to -div(k*grad u)=f of such a polynomial. I have a copy of Understanding The Finite Element Method by Gockenbach’s to refer to but really did not yet find it where such a derivation is explained. It seems this polynomial is manufactured for the unit square described in the comment block. How is the above polynomial @ ( f=vectorize(inline(… ) derived?