For a fixed natural number $n$ is it possible to obtain an asymptotic for the number of solutions $(a,x,b,y)\in\mathbb{N^4}$ to $ax+by=n$ or equivalently an asymptotic for $\sum_{k=1}^{n-1}d(k)d(n-k)$?

Using some heuristics I got that possibly:

$$\sum_{k=1}^{n-1}d(k)d(n-k)\sim \frac{3}{\pi^2}\sigma(n)\ln(n)^2$$

I wrote:

$$\sum_{k=1}^{n-1}d(k)d(n-k)=\sum_{ax+by=n}_{(a,x,b,y)\in \mathbb{N^4}}1=\sum_{d\mid n}\sum_{ax+by=n}_{(a,b)=d}_{(a,x,b,y)\in \mathbb{N^4}}1=\sum_{d\mid n}\sum_{ax+by=\frac{n}{d}}_{(a,b)=1}_{(a,x,b,y)\in \mathbb{N^4}}1$$ $$=\sum_{d\mid n}\sum_{(a,b)\in \mathbb{N^2}}_{(a,b)=1}\sum_{ax+by=\frac{n}{d}}_{(x,y)\in \mathbb{N^2}}1$$

Now I considered a fixed natural number $m$ and looked at the solutions $(x,y)$ to $ax+by=m$

Then noted that $x\leq \frac{m-b}{a}$ and that the solutions $x$ would be of the form: $x\equiv c+a,c+2b,c+3b,c+4b,...c+kb$, so I wrote $c+kb\leq\frac{m-b}{a}\implies k\leq \frac{m-b}{ab}-\frac{c}{b}$

So the largest such $k$ that should satisfy this is $\lfloor{\frac{m-b}{ab}-\frac{c}{a}}\rfloor\approx \lfloor{\frac{m}{ab}}\rfloor$

Then went back and wrote:

$$\sum_{k=1}^{n-1}d(k)d(n-k)\approx \sum_{d\mid n}\sum_{(a,b)\in \mathbb{N^2}}_{(a,b)=1}\lfloor{\frac{n/d}{ab}}\rfloor\approx \sum_{d\mid n}\frac{n}{d}\sum_{ab\leq n}_{(a,b)=1}_{(a,b)\in \mathbb{N^2}}\frac{1}{ab}$$ $$=\sigma(n)\sum_{ab\leq n}_{(a,b)=1}\frac{1}{ab}=\sigma(n)\sum_{k\leq \sqrt{n}}\frac{\mu(k)}{k^2}\sum_{ab\leq \frac{n}{k^2}}\frac{1}{ab}=\sigma(n)\sum_{k\leq \sqrt{n}}\frac{\mu(k)}{k^2}\sum_{j\leq \frac{n}{k^2}}\frac{d(j)}{j}$$ $$\approx\sigma(n)\sum_{k\leq \sqrt{n}}\frac{\mu(k)}{k^2}\frac{\ln(n/k^2)^2}{2}\sim\frac{\sigma(n)\ln(n)^2}{2}\sum_{k\ge 1}\frac{\mu(k)}{k^2}=\frac{3}{\pi^2}\sigma(n)\ln(n)^2$$

But this doesn't seem right either because I know that,

$$\sum_{ax+by\leq m}_{(a,x,b,y)\in \mathbb{N^4}}1=\frac{m^2\ln(m)^2}{2}+(2\gamma-\frac{3}{2})m^2\ln(m)+(2\gamma^2-\frac{5}{2}\gamma-\frac{\pi^2}{12}+\frac{3}{2})m^2+O(m^{1+\theta}\ln(m))$$

Where $\theta$ is the same $\theta$ in Dirichlet's divisor problem.

Thus I would expect that if $\sum_{k=1}^{n-1}d(k)d(n-k)\sim C'\sigma(n)\ln(n)^2$ for some constant $C'$ that we would then have that $C'=\frac{6}{\pi^2}$, which it isn't according to the last heuristic.

So I'm not really sure how to proceed, and would really appreciate any help in finding an asymptotic expansion for $\sum_{k=1}^{n-1}d(k)d(n-k)$.