We consider the functional equation $$\tag1f(x)^2-f(y)^2=f(x+y)f(x-y)$$ for all $x,y\in\mathbb R$.
Using straightforward calculation (esp. the addition theorems) one verifies that the following functions are among the solutions for $(1)$:
- $f(x)=\kappa x$
- $f(x)=\kappa\sin\omega x$
- $f(x)=\kappa\sinh \omega x$
We shall call these solutions Norbert-solutions as Norbert showed in his answer that they make up all $C^2$ solutions of $(1)$.
However, the problem statement was not about twice differentiable solutions, but rather about bounded continuous solutions. We shall see that all bounded continuous solutions are among the Norbert-solutions.
Let $f$ be a bounded continuous solution of $(1)$ that is not identically zero. As has been pointed out $f(0)^2-f(0)^2=f(0)f(0)$ implies that $f(0)=0$.
Let $P\subseteq \mathbb R$ denote the set of periods of $|f|$. Because $f(0)=0$, we have $P\subseteq f^{-1}(0)$. On the other hand, if $a\in f^{-1}(0)$ then $f(a+x)^2-f(x)^2=f(a+2x)f(a)=0$ implies $|f(a+x)|=|f(x)|$ for all $x$, that is $a\in P$. In other words: $P=f^{-1}(0)$. As $P$ is a proper (because $f$ is not identically zero) closed (because $f$ is continuous) subgroup (always) of $\mathbb R$, we have $f^{-1}(0)= p\mathbb Z$ for some $p\ge 0$.
If $a\in P$ then also $2a\in P$, hence
$f(a+x)^2-f(a-x)^2=f(2a)f(2x)=0$ and $f(a+x)f(a-x)=f(a)^2-f(x)^2\le 0 $, i.e. $$ \tag 2f(a-x)=-f(a+x)\qquad\text{if }f(a)=0.$$
If $p>0$, let $I=(0,\frac12p)$; otherwise let $I=(0,\infty)$.
For $h\in I$ consider the sequence $a_k:=f(kh)$, $k\in \mathbb N$. By choice of $I$ and the IVT, $a_1=f(h)$ and $a_2=f(2h)$ have the same sign, i.e. $\frac{a_2}{2a_1}>0$.
The sequence $(a_k)_{k=1}^\infty$ is uniquely determined by $a_1,a_2$ because from $(1)$ and $(2)$ we obtain the recursion
$$\tag3 a_{k+1}=\begin{cases}\frac{a_k^2-a_1^2}{a_{k-1}}&\text{if }a_{k-1}\ne0,\\
-a_{k-3}&\text{if }a_{k-1}=0\end{cases}$$
for $k\ge 2$. (Note that $a_{k-1}= 0$ implies $k\ge 4$).
There exists a Norbert-solution $g\in C^2(\mathbb R)$ with $g(h)=a_1$, and $g(2h)=a_2$ (this uses the "angle duplication formulas" for $\sin$ and $\sinh$; we do not require differential calculus as one might simply verify immediately that these functions are solutions of $(1)$):
- (Parabolic case) If $\frac{a_2}{2a_1}=1$, let $g(x)=\kappa x$ with $\kappa=\frac{a_1}h $.
- (Elliptic case) If $0<\frac{a_2}{2a_1}<1$ let $g(x)=\kappa \sin\omega x$ with $\omega=\frac1h \arccos\frac{a_2}{2a_1}$ and $\kappa=\frac{2a_1\cdot |a_1|}{\sqrt{4a_1^2-a_2^2}}$
- (Hyperbolic case) If $\frac{a_2}{2a_1}>1$ let $g(x)=\kappa \sinh\omega x$ with $\omega=\frac1h \operatorname{arcosh}\frac{a_2}{2a_1}$ and $\kappa=\frac{2a_1\cdot |a_1|}{\sqrt{a_2^2-4a_1^2}}$
By induction, $g(kh)=a_k$ for all $k\in\mathbb N$. In the parabolic and hyperbolic case, the sequence $a_k$ is therefore unbounded, contrary to the assumption that $f$ is bounded.
Thus $a_k=\kappa \sin k\omega h$ for all $k\in\mathbb N$. The angle $\omega h=\arccos\frac{a_2}{2a_1}$ is in $(0\frac\pi 2)$, hence the sine certainly switches signs. Therefore $f$ has nontrivial zeroes, i.e. $p>0$.
A priori, $\kappa$ and $\omega$ depend on $h$. If picking $h'=\frac 12 h$ leads to $\omega', \kappa'$, then $\omega'$ must be "compatible" with $\omega$, that is $\cos h\omega'=\cos h\omega$ and hence $2h'\omega'=h\omega'=\pm h\omega+2m\pi$. By the properties of $\arccos$, we have $0< h'\omega<\frac \pi2$ and $0<h\omega<\frac\pi 2$ so that we can conclude $\omega'=\omega$. Then automatically also $\kappa'=\kappa$.
This means that the Norbert-solution $g(x)=\kappa\sin\omega x$ that we picked for $h$ coincides with $f$ not only at all $x\in h\mathbb Z$, but also for $x\in \frac12h\mathbb Z$ and by induction for all $x\in\frac1{2^n}h\mathbb Z$
and ultimately on the dense set $\bigcup_{n\in\mathbb N}\frac1{2^n}h\mathbb Z$. By continuity, $f=g$.
Remark: We dropped the parabolic and hyperbolic cases only because the problem statement asked for bounded solutions. A quick thought about these cases shows that we again have uniqueness of $\omega$ and $\kappa$ in these cases. Therefore, all continuous solutions of $(1)$ are Norbert-solutions.