Consider a function $f:\mathbb{R}_+^2\rightarrow\mathbb{R}$ of two non-negative real variables (or more generally of several real variables) that is increasing in each argument, continuous, additively (or multiplicatively) separable, that is, it can be written in the form $$f(x,y)=a(x)+b(y)$$ for functions of one variable, $a$ and $b$, and satisfies the following notion of "scale invariance": For each $\lambda>0$ and $(x,y),(\tilde{x},\tilde{y})\in\mathbb{R}_+^2$ we have $$f(x,y)=f(\tilde{x},\tilde{y})\Leftrightarrow f(\lambda x,\lambda y)=f(\lambda \tilde{x},\lambda \tilde{y}).$$ Is there an example of such a function that is not everywhere differentiable? We tried hard to find such an example. For example, we know that if such functions exist, they have to be non-differentiable on a dense null set, by the scale invariance property.
Edit: The original question was asked for non-continuous functions and it was answered negatively. In the edited version we ask the question for continuous functions.
Best Answer
If we allow $f$ to be discontinuous, then the answer is yes: $f$ need not be differentiable.
We can choose $a(x)$ and $b(y)$ so that their ranges $A$ and $B$ are Cantor-like sets which have the following property: $$ \text{if $\alpha, \alpha' \in A$ and $\beta, \beta' \in B$, and $\alpha + \beta = \alpha' + \beta'$, then $\alpha = \alpha'$ and $\beta = \beta'$.} $$ Then $f(x, y) = f(x', y')$ implies $x = x'$ and $y = y'$: just use the above property with $\alpha = a(x)$, $\beta = b(y)$, $\alpha' = a(x')$, $\beta' = b(y')$. This means that the scale-invariance property is trivially satisfied.
Examples of such sets $A$ and $B$ are easy to construct using decimal expansions. For instance, $A$ can be the set of real numbers which can be written using only $0$ and $1$, while $B$ — with $0$ and $2$. Then the function $a(x)$ can be defined as follows: in order to compute $a(x)$, write the binary expansion of $x$ and interpret it as the decimal expansion of $a(x)$. The function $b(y)$ can be defined in a similar way, or simply as $b(y) = 2 a(y)$.
On the other hand, if we require $f$ to be continuous, then the answer is no: $f$ is necessarily differentiable.
Step 1. Fix a positive real $\lambda \geqslant 1$. Since $f(1, 1) \leqslant f(\lambda, 1) \leqslant f(\lambda, \lambda)$, there is a number $\mu = \mu(\lambda) \in [1, \lambda]$ such that $f(\mu, \mu) = f(\lambda, 1)$ (here we use continuity of $f$). By scale-invariance, with $\lambda = x / y$, we have $f(x, y) = f(\lambda y, y) = f(\mu y, \mu y)$. A similar argument works if $0 < \lambda \leqslant 1$. Thus, if we denote $\phi(x) = f(x, x)$, then $$ f(x, y) = \phi(y \mu(\tfrac xy)) . $$
Step 2. If $f(x, y) = a(x) + b(y)$, then the above equality takes form $$ a(x) + b(y) = \phi(y \mu(\tfrac xy)) . $$ In particular, if $x_1 < x_2$, we have $$ a(x_2) - a(x_1) = (a(x_2) + b(x_1)) - (a(x_1) + b(x_1)) = \phi(x_1 \mu(\tfrac{x_2}{x_1})) - \phi(x_1) $$ and $$ a(x_2) - a(x_1) = (a(x_2) + b(x_2)) - (a(x_1) + b(x_2)) = \phi(x_2) - \phi(x_2 \mu(\tfrac{x_1}{x_2})) . $$ If $x_1 = t$ and $x_2 = s t$, then the above identities read $$ a(s t) - a(t) = \phi(t \mu(s)) - \phi(t) = \phi(s t) - \phi(\nu(s) t),$$ where $\nu(\lambda) = \lambda \mu(1/\lambda)$ also lies between $1$ and $\lambda$. Similarly, $$ b(s t) - b(t) = \phi(s t) - \phi(t \mu(s)) = \phi(\nu(s) t) - \phi(t) . $$
Step 3. If $\mu(s) = 1$ for some $s \ne 1$, then $a(s t) - a(t) = 0$, and hence $a$ is constant. Similarly, if $\mu(s) = s$ for some $s \ne 1$, then $\nu(1/s) = 1$, so that $b(t / s) - b(t) = 0$ and hence $b$ is constant. Thus, in what follows we assume that $\mu(s)$ lies strictly between $1$ and $s$ for all $s \ne 1$.
Step 4. We already know that $$ \phi(t \mu(s)) - \phi(t) = \phi(s t) - \phi(\nu(s) t) . $$ If $\psi$ is a nonnegative, smooth, compactly supported function and $$ \Phi(t) = \int_0^\infty \phi(u t) \psi(u) dt , $$ then $\Phi$ is smooth and $$ \Phi(t \mu(s)) - \Phi(t) = \Phi(s t) - \Phi(\nu(s) t) . $$ We consider the Taylor expansion of the above expressions near $s = 1$. Unfortunately, this gets rather technical (I guess a simpler approach is possible, but I fail to see one straight away).
First, we find that $$ \lim_{s \to 1} \frac{s - \nu(s)}{\mu(s) - 1} = \frac{\lim_{s \to 1} \frac{\Phi(t \mu(s)) - \Phi(t)}{\mu(s) - 1}}{\lim_{s \to 1} \frac{\Phi(s t) - \Phi(\nu(s) t)}{s - \nu(s)}} = \frac{t \Phi'(t)}{t \Phi'(t)} = 1 , $$ and similarly $$ \lim_{s \to 1} \frac{s - \mu(s)}{\nu(s) - 1} = 1 . $$ Therefore, $$ \begin{aligned} & \lim_{s \to 1} \frac{(s - 1)^2 - (\nu(s) - 1)^2 - (\mu(s) - 1)^2}{(\mu(s) - 1)(\nu(s) - 1)} \\ & \qquad = \lim_{s \to 1} \frac{(s - \mu(s)) (s - \nu(s))}{(\mu(s) - 1)(\nu(s) - 1)} + \lim_{s \to 1} \frac{(s - \mu(s)) (\mu(s) - 1)}{(\mu(s) - 1)(\nu(s) - 1)} \\ & \qquad \qquad + \lim_{s \to 1} \frac{(s - \nu(s)) (\nu(s) - 1)}{(\mu(s) - 1)(\nu(s) - 1)} - \lim_{s \to 1} \frac{(\mu(s) - 1)(\nu(s) - 1)}{(\mu(s) - 1)(\nu(s) - 1)} \\ & \qquad = 1 + 1 + 1 - 1 = 2 . \end{aligned} $$ It follows that $$ 0 = \lim_{s \to 1} \frac{\Phi(s t) - \Phi(\nu(s) t) - \Phi(t \mu(s)) + \Phi(t)}{(\mu(s) - 1) (\nu(s) - 1)} = p t \Phi'(t) + t^2 \Phi''(t) , $$ where a finite limit $$ p = \lim_{s \to 1} \frac{s - \nu(s) - \mu(s) + 1}{(\mu(s) - 1) (\nu(s) - 1)} $$ necessarily exists. We conclude that $\Phi(t) = c_1 t^{1 - p} + c_2$, unless $p = 1$, in which case $\Phi(t) = c_1 \log t + c_2$.
Since the function $\psi$ was arbitrary, we necessarily have $\phi(t) = C_1 t^{1 - p} + C_2$, unless $p = 1$, in which case $\Phi(t) = C_1 \log t + C_2$.
Step 5. Recall that $$ a(s t) - a(t) = \phi(t \mu(s)) - \phi(t) .$$ Dividing both sides by $(\mu(s) - 1)$ and passing to the limit as $s \to 1$, we find that $$ \lim_{s \to 1} \frac{a(s t) - a(t)}{\mu(s) - 1} = t \phi'(t) . $$ Since this holds for every $t$, the function $a$ is necessarily differentiable, with $$ \frac{t a'(t)}{\mu'(0)} = t \phi'(t) . $$ Thus, $a(t) = C_3 + \mu'(0) \phi(t)$. Similarly, $b(t) = C_4 + \nu'(0) \phi(t)$.
Summary. We have shown that $a(x)$ and $b(y)$ are differentiable, and in fact for some constants $c_1, c_2, c_3, c_4, \gamma$, with $c_2 \gamma, c_4 \gamma > 0$ if $\gamma \ne 0$, we have $$ a(x) = c_1 + c_2 x^\gamma, \qquad b(y) = c_3 + c_4 y^\gamma , $$ except for the case $\gamma = 0$, where we have $c_2, c_4 > 0$ and $$ a(x) = c_1 + c_2 \log x, \qquad b(y) = c_3 + c_4 \log y . $$