$\renewcommand{\Re}{\operatorname{Re}}$ $\renewcommand{\Im}{\operatorname{Im}}$ $\newcommand{\bR}{\mathbb{R}}$ $\newcommand{\bC}{\mathbb{C}}$ $\newcommand{\bZ}{\mathbb{Z}}$ $\newcommand{\const}{\operatorname{const}}$
We are interested in solving problem \begin{align} &P(x,-ih \nabla , h)u=0, \label{eq-5.4.1}\\ &(-ih\partial_t)^j u|_{t=0}= f_j(x')&& j=0,\ldots, m-1 \label{eq-5.4.2} \end{align} where $x=(x_0,x_1,\ldots, x_d)$, $t=x_0$, $x'=(x_1,\ldots, x_d)$, $m$ is maximal degree of $-ih\partial_t$ in $P$ (and also in $P_0$); other derivatives can have larger degrees.
For Schrödinger equation (5.1.3) $m=1$, for wave equation (5.1.1) $m=2$.
We assume that initial functions $f_k(x')$ are of the form \begin{equation} f_j(x')\sim e^{ih^{-1}S_0(x')}\sum_{n\ge 0} F_{j,n} (x') h^n \label{eq-5.4.3} \end{equation} (the same phase but different amplitudes). We are looking for solution in the form \begin{equation} u(x)\sim \sum_{1\le k\le m} e^{ih^{-1}S_k(x)}\sum_{n\ge 0}A_{n,k} (x) h^n. \label{eq-5.4.4} \end{equation} We know that $S_k(x)$ must satisfy (5.2.1) \begin{equation} P_0(x,\nabla S_k(x))=0. \label{eq-5.4.5} \end{equation}
Condition 1. All roots $p_0=-H_k (x,p')$, $k=1,\ldots, m$ of the polynomial $P_0(x,p_0, p')$ are real and distinct as $p'=\nabla' S_0(x')$.
Obviously this condition is fulfilled for Schrödinger equation, and it is fulfilled for wave equation as $\nabla' S_0(x')\ne 0$.
Under this condition we can find (locally) functions $S_k(x)$ satisfying Hamilton-Jacobi (5.2.3) \begin{equation} \partial_tS_k + H_k(x,\nabla 'S_k)=0\qquad k=1,\ldots,m. \label{eq-5.4.6} \end{equation} Then amplitudes $A_{k,n}$ must satisfy transport equations (5.3.1) (as $n=0$) and similar \begin{equation} \bigl(\partial_t+\sum_{1\le j\le d} H_k^{(j)} (x, \nabla S(x))\partial_j +i Q_k(x)\bigr)A_{k,n} (x)= \sum_{l\le n-1} \mathcal{L}_{k,n-l} A_{k,l}. \label{eq-5.4.7} \end{equation} Here $A_{k,n}$ are defined recurrently. To define $A_{k,n}$ we need their values as $t=0$ (and then we use transport equations).
Consider first $n=0$. Then (\ref{eq-5.4.2})--(\ref{eq-5.4.3}) are fulfilled modulo $O(h)$ if \begin{equation} \sum_{1\le k\le m} (\partial_t S_k)^j A_{k,0}\bigl|_{t=0}=F_{j,0}\qquad j=0,\ldots, m-1. \label{eq-5.4.8} \end{equation} This is a linear $m\times m$ system with Vandermonde matrix and the determinant \begin{equation} \prod _{1\le j< k\le m} \bigl(-H_k (x',\nabla'S_0)+ H_k (x',\nabla'S_0)\bigr)\ne 0 \label{eq-5.4.9} \end{equation} where we used (\ref{eq-5.4.6}) and Condition 1.
Assume that we defined $A_{k,l}$ for all $l=0,\ldots, n-1$ and all $k=1,\ldots, m$ so that initial conditions are fulfilled modulo $O(h^n)$. Then for $A_{k,n}$ we get \begin{equation} \sum_{1\le k\le m} (\partial_t S_k)^j A_{k,n}\bigl|_{t=0}=G_{j,n}\qquad j=0,\ldots, m-1 \label{eq-5.4.10} \end{equation} where $G_{j,n}$ are $F_{j,n}$ and what came from $A_{k,l}$ with $l=0,\ldots, n-1$, $k=1,\ldots, m$.
Again ddeterminant does not vanish.
Consider specifically wave equation (5.1.1) in domain $X$: \begin{equation} u_{tt}-\nabla \cdot (c^2(x)\nabla u)=0 \label{eq-5.4.11} \end{equation} There is a solution \begin{equation} u\sim e^{ih^{-1}\phi (x)}\sum_{n\ge 0} A_k (x)h^n \label{eq-5.4.12} \end{equation} where we prefer to denote eikonal by $\phi$ (etc). But this solution does not satisfy Dirichlet boundary condition \begin{equation} u|_Y=0 \label{eq-5.4.13} \end{equation} or Neumann or Robin boundary condition \begin{equation} (\boldsymbol{\nu}\cdot \nabla +\kappa) u|_Y=0 \label{eq-5.4.14} \end{equation} where $Y=\partial X$ and $\boldsymbol{\nu}$ is a normal to $Y$ directed into $X$.
Assumption 2. It is incoming wave which means that as $t$ increases trajectories come from $X$ to $Y$. This is equivalent to \begin{equation} \frac{\partial \phi}{\partial\boldsymbol{\nu}} : \frac{\partial \phi}{\partial t} \Bigl|_Y >0. \label{eq-5.4.15} \end{equation} where $\frac{\partial \phi}{\partial\boldsymbol{\nu}} =\boldsymbol{\nu}\cdot \nabla$.
To fulfill boundary condition we add a reflected wave \begin{equation} u\sim e^{ih^{-1}\phi (x)}\sum_{n\ge 0} A_k (x)h^n + e^{ih^{-1}\psi (x)}\sum_{n\ge 0} B_k (x)h^n. \label{eq-5.4.16} \end{equation} Both $\phi$ and $\psi$ must satisfy eikonal equation \begin{equation} \phi_t^2 = c^2|\nabla \phi|^2,\qquad \psi_t^2 = c^2|\nabla \psi|^2; \label{eq-5.4.17} \end{equation} they must coincide of $Y$: \begin{equation} \psi=\phi \qquad \text{on }\ Y \label{eq-5.4.18} \end{equation} and differ. Then \begin{equation} \frac{\partial \psi}{\partial\boldsymbol{\nu}}= -\frac{\partial \phi}{\partial\boldsymbol{\nu}} \text{on }\ Y \label{eq-5.4.19} \end{equation} Then the added term is an outgoing wave \begin{equation} \frac{\partial \psi}{\partial\boldsymbol{\nu}} : \frac{\partial \psi}{\partial t} \Bigl|_Y < 0. \label{eq-5.4.20} \end{equation} From this and equations for rays follows the well known rule: reflection angle equals incidence angle.
Remark 1. One can understand this from toy-model $c=\const$ and $X=\{x_1>0\}$, $\phi = ct - k x_1- l x_2$; then $\psi = ct + k x_1- l x_2$ ($k^2+ l^2=1$).
Reflection: $\alpha_{\text{inc}}=\alpha_{\text{refl}}$
Now Dirichlet or Robin boundary condition imply that \begin{equation} B_0= \mp A_0\quad\text{on }\ Y \label{eq-5.4.21} \end{equation} and similar conditions for $B_n$. Since $B_n$ must satisfy transport equations we can find them locally.
Consider now two wave equations (5.1.1) in domains $X_1$ and $X_2$: \begin{equation} u_{j,tt}-\nabla \cdot (c_j^2(x)\nabla u_j)=0 \quad\text{in }\ X_j. \label{eq-5.4.22} \end{equation} These domains have a common boundary $Y$ where two boundary conditions connecting $u_1$ and $u_2$ must be satisfied. We consider incoming wave (\ref{eq-5.4.12}) in $X_1$ satisfying (\ref{eq-5.4.15}) where $\boldsymbol{\nu}$ is a normal to $Y$ directed into $X_1$.
To satisfy two boundary conditions on $Y$ we need not only consider reflected wave (so solution is given by (\ref{eq-5.4.16}) in $X_1$ but also a refracted wave \begin{equation} u\sim e^{ih^{-1}\chi (x)}\sum_{n\ge 0} C_k (x)h^n\quad\text{in }\ X_2. \label{eq-5.4.23} \end{equation} Here $\chi$ must satisfy eikonal equation $\chi_t^2=c_2^2|\nabla \chi|^2$ and correspond to outgoing wave. From this and equations for rays follows the well known Snell law: \begin{equation} \frac{\sin(\alpha_1)}{c_1}=\frac{\sin(\alpha_2)}{c_2} \label{eq-5.4.24} \end{equation} where $\alpha_1 $ is an angle between incidental angle and normal and $\alpha_2$ is an angle between reflection angle and normal.
Remark 2. One can understand this from toy-model $c_j=\const$ and $X_1=\{x_1>0\}$, $X_2=\{x_1< 0\}$, $\phi = c_1t - k x_1- l x_2$; then $\psi = c_2 + k x_1- l x_2$ ($k^2+ l^2=1$) but $\chi = c_1 t +m x_1 -l x_2$; then $k^2+l^2=1$ but $m^2 + l^2= c_1^2/c_2^2$. Clearly $l=\sin (\alpha_1)=\sin (\alpha_2) c_2/c_1$.
Reflection and refraction: $\alpha_{\text{inc}}=\alpha_{\text{refl}}$, $\dfrac{\sin(\alpha_{\text{inc}})}{c_1}=\dfrac{\sin(\alpha_{\text{refr}})}{c_2}$
Remark 3. It may happen that we cannot find real $\alpha_2$ as \begin{equation} \frac{c_2\sin(\alpha_1)}{c_1} >1. \label{eq-5.4.25} \end{equation} Then we cannot find real valued $\chi$ but we can find complex-values $\chi$ with $\Im \chi > 0$. Then we have a wave which exponentially decays into $X_2$ penetrating there on the depth $\asymp h$. Geometrically this wave is not observable (it is where incoming wave hits $Y$) but analytically it is still here.
This is called a complete internal reflection.
Remark 4. Propagation theory near boundary near tangency of trajectories to the boundary is much more complicated and not completely explored. Even two easiest cases are really difficult:
Remark 5. Another very difficul case is scattering on polygons and polyhedra or cones
Pictures will be here, or may be on the lectures only
Remark 6.