I am not at all an expert in this field. However, I believe I can show that there are no examples where $Df$ vanishes to finite order. Thus, in particular there are no analytic solutions and, if the OP is right about two-jet-determination, there are no solutions.
Let $f: M \to N$ be a weakly conformal map between Riemannian manifolds of dimension $d>2$. We write $g_M$ and $g_N$ for the metrics on $M$ and $N$, so the conformal condition is that $(Df)^T g^N (Df) = \lambda g^M$ for some scalar valued function $\lambda$ on $M$. If $f$ is singular at $x_0 \in M$, then we must have $\lambda(x_0)=0$, so we must have $Df=0$. We say that $f$ vanishes to order $k$ at $x_0$ if all $r$-fold derivatives of $f$ vanish at $x_0$ for $r<k$ from some $k$-fold derivative doesn't. So we have shown that, if $f$ is singular at $x_0$, then it vanishes to some order $k>1$. We are going to show that there are no solutions which vanish to finite order $>1$.
The problem is local, so we can assume that $M$ and $N$ are open subsets of $\mathbb{R}^d$ and we may take $x_0=0$ and $f(x_0)=0$. Furthermore, we may take bases such that $g^M(0)$ and $g^N(0)$ are the standard inner product. These choices still allow us the ability to independently make further orthogonal changes of basis on source and target. We expand $f$ in a Taylor series as $f(x) = f_k(x) + O(|x|^{k+1})$, expand the metrics as $g^M(x) = \mathrm{Id}+ O(|x|)$ and $g^N(y) = \mathrm{Id}+O(|y|)$, and expand the scalar $\lambda$ as $\lambda(x)=\lambda_m(x) +O(|x|^{m+1})$. Here $f_k$ and $\lambda_m$ are nonzero polynomials of degree $k$ and $m$.
The conformal condition is that $$(Df_k + \cdots)^T (\mathrm{Id}+ \cdots) (Df_k + \cdots) = (\lambda_m + \cdots ) (\mathrm{Id}+ \cdots).$$ Since $f_k \neq 0$, the matrix $Df_k$ is nonzero and thus $(Df_k)^T (Df_k)$ is nonzero. So the leading terms $(Df_k)^T (Df_k)$ and $\lambda_m \mathrm{Id}$ must have matching degrees. We deduce that $m = 2k-2$ and $(Df_k)^T (Df_k) = \lambda_{m} \mathrm{Id}$.
Thus, the polynomial map $f_k$ is conformal between $\mathbb{R}^d$ and itself, with the standard Euclidean structure. This contradicts the result of Iwaniec and Martin cited by the OP.
Also, in the case $k=2$, we can verify that there is no solution by a direct computation. I will rename my variables -- my goal is to show that there is no conformal map $f: \mathbb{R}^d \to \mathbb{R}^d$ given by quadratic polynomials such that $(Df)(x)^T (Df)(x) = \lambda(x) \mathrm{Id}$ where $\lambda$ is a polynomial of degree $2$. Let the components of $f$ be $f(x) = \tfrac{1}{2} (x^T P_1 x, x^T P_2 x, \cdots, x^T P_n x)$ and let $\lambda(x) = \tfrac{1}{2} x^T Q x$, where the $P_i$ and $Q$ are symmetric.
We can break the matrix $Df$ up into rows: $$Df = \begin{bmatrix} x^T P_1 \\ x^T P_2 \\ \vdots \\ x^T P_n \end{bmatrix}.$$ Our given condition is that $(Df)(x)^T (Df)(x) = \lambda(x) \mathrm{Id}$; since $\lambda(x)$ is nonzero, we can work in the field of rational functions of $x$ and deduce that we also have $(Df)(x) (Df)(x)^T = \lambda(x) \mathrm{Id}$, giving $$\begin{bmatrix} x^T P_1 \\ x^T P_2 \\ \vdots \\ x^T P_n \end{bmatrix} \begin{bmatrix} P_1 x & P_2 x & \cdots & P_n x \end{bmatrix} = \tfrac{1}{2} (x^T Q x) \mathrm{Id}.$$ So, for all $x \in \mathbb{R}^d$, we have $$x^T P_i P_j x = \begin{cases} \tfrac{1}{2} x^T Q x & i = j \\ 0 & i \neq j \end{cases}.$$ We can't deduce that $P_i P_j = \tfrac{1}{2} Q$ from this, since $P_i P_j$ may not be symmetric. But we do deduce $$P_i P_j + P_j P_i = \begin{cases} Q & i=j \\ 0 & i \neq j \end{cases}. \quad (\ast)$$
Making an orthogonal change of basis, we may assume that $Q$ is diagonal. Let the eigenvalues of $Q$ be $q_1$, $q_2$, ... $q_k$ with multiplicities $d_1$, $d_2$, ..., $d_k$. Since $Q = 2 P_i^2$, all the $P_i$ commute with $Q$ and thus they are block diagonal with blocks of size $d_1$, $d_2$, ..., $d_k$. Thus, we may solve the equation $(\ast)$ separately in each block, and may thus assume that $Q$ is diagonal; say $Q = q \mathrm{Id}$. Since $Q = 2 P_i^2$, we have $q>0$; rescaling each $P$ by $\tfrac{1}{\sqrt{2 q}}$, we may assume that $Q = \mathrm{Id}$.
Thus, $(\ast)$ becomes $$P_i P_j + P_j P_i = \begin{cases} \mathrm{Id} & i=j \\ 0 & i \neq j \end{cases}. \quad (\ast)$$ In other words, $\mathbb{R}^{d_i}$ is a representation of the Clifford algebra on $d$ generators with $d_i \leq d$. For $d>2$, this is impossible.