We can explicitly give the requested implicit definition for multiplication: It is the unique function on $(\mathbb{N},S,0)$ satisfying:
\begin{align}
0a&=0\\
(S0)a&=a\\
ab&=ba\\
a(bc)&=(ab)c\\
(Sa)b=(Sa)c &\to b=c\\
(Sa)Sa &= S(aSSa)\\
(SS0)Sa &= SS(aSS0)\\
(S(aSS0))(S(ab))&= S(aSS(bS(aSS0)))
\end{align}
The last three identities would be more familiarly written as
\begin{align}
(a+1)^2 &= 1+a(2+a)\\
2(1+a)&=2+2a\\
(1+2a)(1+ab)&=1+a(2+b(1+2a))
\end{align}
We refer to them as the *square*, *even* and *odd* identities: they give the key facts about squares and multiplication by even and odd numbers. The square identity defines the square of $a+1$ in terms of multiplication by $a$, the odd identity defines (a special case of) multiplication by $2a+1$ in terms of multiplication by $a$, and we can use them in inductions. 

(The odd identity is a special case of $1+a(1+b(1+c(1+a(1+b))))=(1+a+ab)(1+abc)$, which may be easier to understand.)

As is usual in these matters, we look at numerals of the form $S^n0$. For each positive integer $n$, that numeral is a term in the language of the model. We quantify over $n$'s outside the model.

We prove by strong induction that for all $m$ and $n$, the axioms imply $S^m0\ S^n0=S^{mn}0$. This determines the values of the multiplication function for the whole domain of the model.

- The cases $m=0$ and $m=1$ are immediate.

- The case $m=2$ follows as an induction from the even identity.

- The composite case $m=jk$ with $j<m$, $k<m$ is $$S^{jk}0\ S^n0 = S^j0\ S^k0\ S^n0 = S^j0\ S^{kn}0 = S^{jkn}0$$ which follows from the cases $m=j$, $m=k$ and $m=j$. 

- This leaves the case where $m$ is an odd prime, and we take $m=1+2j$. 

- - If $n=0,1,2$ or $n=jk$ then we use the same argument as for $m$, so we may assume wlog that $n$ is also an odd prime.

- - If $m=n$, then
$$S^m0\ S^m0 = S(S^{m-1}0\ S^{m+1}0)=S(S^{m^2-1}0)=S^{m^2}0$$
by the square identity, the inductive hypothesis for $m-1$, and simplifying notation. So we may also assume wlog that $n<m$.

- - If $n=1+jk$ (the mod-1 case), then abbreviating $a=S^j0$, $b=S^k0$ gives
\begin{align}
S^{1+2j}0\ S^{1+jk}0
&=S(aSS0)\ S(ab) \ \ \ \text{by inductive hypotheses}\\
&=S(aSS(bS(aSS0))) \ \ \ \ \text{by the odd identity}\\
&=S(aSS(bS(S^{2j}0))) \ \ \ \ \text{by inductive hypothesis for }m=j\\
&=S(aSS(b(S^{1+2j}0))) \ \ \ \ \text{by simplifying notation}\\
&=S(aSS(S^{(1+2j)k}0)) \ \ \ \ \text{by inductive hypothesis for }m=k\\
&=S(aS^{2+(1+2j)k}0) \ \ \ \ \text{by simplifying notation}\\
&=S(S^{2j+(1+2j)jk}0) \ \ \ \ \text{by inductive hypothesis for }m=j\\
&=S^{(1+2j)(1+jk)}0 \ \ \ \ \text{by simplifying notation}\\
\end{align}

- - If $n=i+jk$, where $1<i<j$, then because $n$ is prime, $i$ must be relatively prime to $j$. So let $h$ be a multiplicative inverse of $i$ mod $j$, with $hi=gj+1$. Then
\begin{align}
S^{1+2j}0\ S^{h(i+jk)}0 &= S^{(1+2j)h(i+jk)}0\ \ \ \ \text{by the mod-1 case}\\ 
S^{1+2j}0\ S^{i+jk}0 &= S^{(1+2j)(i+jk)}0 \ \ \ \ \text{by cancellation}\\
S^m0\ S^n0 &= S^{mn}0
\end{align}

- - In any case $S^m0\ S^n0 = S^{mn}0$, QED.