approximating a maximum function by a differentiable function

Yes it is. One possibility is the following: Note that $\def\abs#1{\left|#1\right|}$ $$ \max\{x,y\} = \frac 12 \bigl( x+ y + \abs{x-y}\bigr), $$ take a differentiable approximation of $\abs\cdot$, for example $\def\abe{\mathop{\rm abs}\nolimits_\epsilon}$$\abe \colon \mathbb R \to \mathbb R$ for $\epsilon > 0$ given by $$ \abe(x) := \sqrt{x^2 + \epsilon}, \quad x \in \mathbb R $$ and define $\max_\epsilon \colon \mathbb R^2 \to \mathbb R$ by $$ \max\nolimits_\epsilon(x,y) := \frac 12 \bigl( x+y+\abe(x-y)\bigr). $$ Another possibility is to take a smooth mollifier $\phi_\epsilon$ and let $\max'_\epsilon :=\mathord\max * \phi_\epsilon$.

Another possibility is given by: \begin{equation} \max (x,y) \approx \ln(e^x+e^y) \end{equation} This approximation doesn't work well for similar values of $x$ and $y$. We can however, remedy this by introducing a scaling parameter $N$: \begin{equation} \max (x,y) \approx \frac1N\ln(e^{N x}+e^{N y}) \end{equation} for large values of $N$.

A general definition is given by:

\begin{equation} \max\limits_{x \in S} \approx \frac1N\ln(\sum\limits_{x\in S} e^{N x}) \end{equation}

Note that in practice $e^{Nx}$ will give unworkably large answers for large $N$.