For 2nd year engineering students, you should not be wasting time trying to prove the expansion, except perhaps giving a heuristic argument for the quadratic approximation in the case of two variables. Stick to operational algebraic formulations and computations, which is the only thing they'll ever need.
My go-to book for an elementary treatment of this topic (indeed, for an elementary treatment of many multivariable calculus topics) is Calculus of Several Variable by Serge Lang (see the section titled "The general expression for Taylor's formula"). Don't let Lang's reputation in his many other books mislead you --- this book is excellent for students with a 2-semester background in single variable calculus. The examples are excellent, the explanatory material is well-written and comprehensible by students, the text is uncluttered both verbally and symbolically, and while the text is suitable for average level students who do not intend to major in math, future math majors will not be misled by incorrect or misleading statements.
What follows is taken from a take-home assignment I gave to very strong high school students here in a post-BC calculus class that I taught several times in the mid to late 1990s. I didn't cover this topic in class (I "worked it in" by including it in a take-home assignment), but I can easily imagine using it as the basis for one class session devoted to this topic, and thus it might be of use to you. Of course, some simple examples would need to be included and probably some aspects of the problems below would not be appropriate in your case. Keep in mind these were near olympiad-level students, so I tried to work in enrichment topics as much as I could, something you should NOT do unless the enrichment topics are engineering applications.
One formulation of the two-variable Taylor expansion is:
$$ \begin{array} ff(x,y) & = & \; & f(a,b) \\ \; & \; & + & \left[\,f_x(a,b) \cdot (x-a) \; + \; f_y(a,b) \cdot (y-b) \right] \, \\ \; & \; & + & \frac{1}{2!} \cdot \left[\,f_{xx}(a,b) \cdot (x-a)^2 \; + \; 2f_{xy}(a,b) \cdot (x-a)(y-b) \; + \; f_{yy}(a,b) \cdot (y-b)^2 \right] \, \\ \; & \; & + & \cdot \cdot \cdot
\end{array} $$
This corresponds to the following formulation of the one-variable Taylor expansion:
$$ f(x) \;\; = \;\; f(a) \;\; + \;\; f'(a) \cdot (x-a) \;\; + \;\; \frac{1}{2!} \cdot f''(a) \cdot (x-a)^2 \;\; + \;\; \cdot \cdot \cdot $$
In this formulation of the one-variable Taylor expansion, if we replace $(x-a)$ with $\Delta x$ (hence, $x$ is replaced with $a + \Delta x),$ then we get an alternative formulation of the one-variable Taylor expression:
$$ f(a + \Delta x) \;\; = \;\; f(a) \;\; + \;\; f'(a) \cdot \Delta x \;\; + \;\; \frac{1}{2!} \cdot f''(a) \cdot \left(\Delta x \right)^2 \;\; + \;\; \cdot \cdot \cdot $$
The corresponding alternative formulation of the two-variable Taylor expansion is:
$$ \begin{array} ff\left(a + \Delta x, \; b + \Delta y\right) & = & \; & f(a,b) \\ \; & \; & + & \left[\,f_x(a,b) \cdot \Delta x \; + \; f_y(a,b) \cdot \Delta y \right] \, \\ \; & \; & + & \frac{1}{2!} \cdot \left[\,f_{xx}(a,b) \cdot \left(\Delta x \right)^2 \; + \; 2f_{xy}(a,b) \cdot \left(\Delta x \right) \left( \Delta y \right) \; + \; f_{yy}(a,b) \cdot \left( \Delta y \right)^2 \right] \, \\ \; & \; & + & \cdot \cdot \cdot
\end{array} $$
Using vector notation, we can express the alternative formulation of the two-variable Taylor expansion in a way that more closely parallels the alternative formulation of the one-variable Taylor expansion. Let
$$ \vec{X_{\,}} = \pmatrix{x \\ y}, \;\;\;\; \vec{X_0} = \pmatrix{a \\ b}, \;\;\;\; \Delta \vec{X_{\,}} = \pmatrix{\Delta x \\ \Delta y} $$
Then
$$ \begin{array} f\left[f\right]_{\vec{X_0} \; + \; \Delta \vec{X_{\,}}} & = & \; & \left[f\right]_{\vec{X_0}} \\
\; & \; & + & \left[ \left( \vec{\nabla \;} \cdot \Delta \vec{X_{\,}} \right) (f) \right]_{\vec{X_0}} \, \\ \; & \; & + & \frac{1}{2!} \left[ \left( \vec{\nabla \;} \cdot \Delta \vec{X_{\,}} \right)^2 (f) \right]_{\vec{X_0}} \, \\
\; & \; & + & \frac{1}{3!} \left[ \left( \vec{\nabla \;} \cdot \Delta \vec{X_{\,}} \right)^3 (f) \right]_{\vec{X_0}} \, \\ \; & \; & + & \cdot \cdot \cdot \, ,
\end{array} $$
where the subscripts refer to evaluations and $ \vec{\nabla \;} \cdot \Delta \vec{X_{\,}} $ represents the (formal) dot product of $\vec{\nabla \;}$ with $\Delta \vec{X_{\;}}.$ The various powers of $ \vec{\nabla \;} \cdot \Delta \vec{X_{\,}} $ are to be expanded "algebraically", meaning (for example) that the square of $\frac{\partial}{\partial x}$ represents $\frac{{\partial}^2}{\partial x^2},$ and so on in the obvious manner. I say "obvious", since you already have the expansion written out in terminology you should understand. Just make the obvious identifications needed for the expansions to be equal. The nice thing about writing the Taylor expansion this way is that if you decide to make everything three variables (or more), then the vector formulation remains the same.
Problem 1
Let $f$ be a real-valued function of the three variables $x,$ $y,$ and $z,$ and put
$$ \vec{X_{\,}} = \pmatrix{x \\ y \\ z}, \;\;\;\; \vec{X_0} = \pmatrix{a \\ b \\ c}, \;\;\;\; \Delta \vec{X_{\,}} = \pmatrix{\Delta x \\ \Delta y \\ \Delta z} $$
Expand the following to obtain three-variable terms that correspond to the two-variable terms in the first two-variable Taylor expansion I gave.
$$ \left[ \left( \vec{\nabla \;} \cdot \Delta \vec{X_{\,}} \right) (f) \right]_{\vec{X_0}} \;\; = \;\; \text{?} $$
$$ \left[ \left( \vec{\nabla \;} \cdot \Delta \vec{X_{\,}} \right)^2 (f) \right]_{\vec{X_0}} \;\; = \;\; \text{?} $$
$$ \left[ \left( \vec{\nabla \;} \cdot \Delta \vec{X_{\,}} \right)^3 (f) \right]_{\vec{X_0}} \;\; = \;\; \text{?} $$
For the last expression, you will have to cube a trinomial. Here is a neat way to obtain the cube of a trinomial.
$(A + B + C)^3 = (A+B+C)(A+B+C)(A+B+C)$ can be algebraically expanded by adding together all products of the form $PQR,$ where $P$ is a choice of a term in the left factor, $Q$ is a choice of a term in the middle factor, and $R$ is a choice of a term in the right factor. Thus, there will be $(3)(3)(3) = 27$ terms added together, many of which will be the same (e.g. choosing $P = A,$ $Q = A,$ $R = B$ gives the term $A^2B$; choosing $P = A,$ $Q = B,$ $R = A$ also gives the term $A^2B).$ As a consequence of this method of expansion, note that all the terms in the algebraic expansion will be of the third degree. These $27$ terms can be organized into the following three types: terms that are the cube of a variable, terms that include the square of a variable, and the terms equal to $ABC.$ Clearly, there is only one term equal to $A^3,$ and similarly for $B^3$ and $C^3.$ Moreover, it is easy to see that there is a total of $3$ terms equal to $A^2B,$ since there is a total of $3$ ways to choose exactly one $B$ and exactly two $A$'s. Similarly for the terms equal to $AB^2,$ $AC^2,$ $A^2C,$ etc. Therefore,
$$ (A+B+C)^3 \; = \; A^3 + B^3 + C^3 + 3(A^2B + A^2C + AB^2 + B^2C + AC^2 + BC^2) + k \cdot ABC $$
for some positive integer constant $k.$ You can find the value of $k$ by substituting $A=B=C=1$ and solving for $k.$
Problem 2
Use the 2nd order Taylor polynomial for functions of $3$ variables to find the quadratic approximation to
$$ f(x,y,z) \; = \; e^x \sin y \cos z \;\;\; \text{about the point} \;\;\; \pmatrix{a \\ b \\ c} = \pmatrix{0 \\ 0 \\ 0} . $$
Problem 3
Use the angle addition formula for SINE to show that $\; \sin y \cos z = \frac{1}{2}\sin(y-z) \; + \; \frac{1}{2} \sin(y+z).$ This is one of the product to sum formulas in trigonometry, which you might recall using in carrying out certain integrations. Now use this equality, and substitution into the one-variable Taylor series for the SINE function, to obtain an expansion of $\sin y \cos z$ (up to and including third degree terms will be enough for what follows). Then multiply this expansion of $\sin y \cos z$ by the 2nd order Taylor polynomial of $e^x,$ namely $1 + x + \frac{1}{2}x^2,$ keeping only terms of degree $2$ or less. Show that the answer you get is the same as what you got in Problem 2.
