Well, not much can be said about the solutions in the general case where $ \kappa $ is arbitrary. The only useful (and fairly trivial) observation is that using induction, you can show that for each $ n \ge 1 $ that is a power of $ 2 $, you have
$$ \delta ( n ) = \delta ( 1 ) - \sum _ { k = 0 } ^ { \log _ 2 ( n ) - 1 } \kappa \left( 2 ^ k \right) \text . \tag 0 \label 0 $$
This shows that given $ \delta ( 1 ) $, $ \delta ( n ) $ is determined uniquely. It very much depends on what $ \kappa $ is whether there is a closed form for $ \delta $ or not.
In the special case that
$$ \kappa ( n ) = \frac M { 8 ( 2 n - 1 ) } \tag 1 \label 1 $$
which is of your interest, WolframAlpha gives the solution in terms of the $ q $-digamma function $ \psi _ q $, which is the $ q $-analog of the digamma function $ \psi $. $ \psi _ q ( x ) $ is defined to be equal to $ \frac 1 { \Gamma _ q ( x ) } \frac { \partial \Gamma _ q ( x ) } { \partial x } $, as we similarly had $ \psi ( x ) = \frac 1 { \Gamma ( x ) } \frac { \mathrm d \Gamma ( x ) } { \mathrm d x } $, where $ \Gamma $ is the gamma function, and $ \Gamma _ q $ is its $ q $-analog, the $ q $-gamma function. While these are all well-know functions, their definitions involve sums, and calculating their values may not be much easier than calculating the sum appearing in \eqref{0}. But as they are well-known special functions, any usual computational software has some implementation of them, which may be useful to you.
Anyways, in case of \eqref{1}, the solution given by WolframAlpha can be simplified to the following form:
$$ \delta ( n ) = K + \frac M 8 \Big( \log _ 2 ( n ) - \log _ 2 ( e ) \cdot \psi _ 2 \big( 1 + \log _ 2 ( n ) \big) \Big) \text , $$
where $ e $ is Napier's constant and $ K $ is a constant which can be determined from the value of $ \delta ( 1 ) $.
In fact, it's not true that linear functions are the only solutions. For any $ a , b \in \mathbb R _ { 0 + } $, the function $ f : \mathbb R \to \mathbb R $ defined with
$$ f ( x ) =
\begin {cases}
a x & x \ge 0 \\
b x & x < 0
\end {cases} $$
is another solution. This can be verified easily, using the fact that $ y $, $ a y $, $ b y $, $ \left( x ^ 2 + a \right) y $ and $ \left( x ^ 2 + b \right) y $ will always have the same sign (in the nonstrict sense), because $ a , b \ge 0 $. Considering these solutions together with linear functions, one can show that there is no other solution, which we will prove.
We continue the argument put forward by the OP, except that we only prove $ g ( y ) = 0 $ for $ y > 0 $, as it may not be true for $ x < 0 $ in the case of the above solutions. Fixing $ y > 0 $, we can choose $ x $ large enough, namely $ x \ge \sqrt { \frac { | a - a y - g ( y ) | } y } $, so that we have $ \left( x ^ 2 + a \right) y + g ( y ) \ge a $, and hence $ g \Big( \left( x ^ 2 + a \right) y + g ( y ) \Big) = 0 $. Then we can use the functional equation for $ g $ to get $ x ^ 2 g ( y ) = - g \big( a y + g ( y ) \big) $, which shows that $ g ( y ) = 0 $, as this is true for all large enough $ x $.
For the negative points, we can take a similar path. Letting $ b = - f ( - 1 ) $ and defining $ h : \mathbb R \to \mathbb R $ with $ h ( x ) = f ( x ) - b x $ for all $ x \in \mathbb R $, we have $ h ( 0 ) = h ( - 1 ) = 0 $ and the functional equation
$$ h \Big( \left( x ^ 2 + b \right) y + h ( y ) \Big) = x ^ 2 h ( y ) + h \big( b y + h ( y ) \big) $$
for all $ x , y \in \mathbb R $. Setting $ y = - 1 $, we get $ h ( x ) = h ( - b ) $ for all $ x \le - b $. In case $ b \le 1 $, we have $ - 1 \le - b $ and since $ h ( - 1 ) = 0 $, $ h ( - b ) = 0 $. In case $ b > 0 $, letting $ y = - b $ and choosing $ x $ large enough so that $ \left( x ^ 2 + b \right) y + h ( y ) \le - b $, we get $ x ^ 2 h ( - b ) = h ( - b ) - h \big( - b ^ 2 + h ( - b ) \big) $, which gives $ h ( - b ) = 0 $ since it's true for all large enough $ x $. So in any case, we have $ h ( x ) = 0 $ for all $ x \le - b $. Again, for any $ y < 0 $, choosing $ x $ large enough so that $ \left( x ^ 2 + b \right) y + h ( y ) \le - b $, we get $ x ^ 2 h ( y ) = - h \big( b y + h ( y ) \big) $, which shows that $ h ( y ) = 0 $, as that is true for all large enough $ x $.
So we know that $ g ( x ) = 0 $ for all $ x \ge \min ( a , 0 ) $ and that $ h ( x ) = 0 $ for all $ x \le \max ( - b , 0 ) $. Note that by definition of $ g $ and $ h $ we have $ g ( x ) - h ( x ) = ( a - b ) x $ for all $ x \in \mathbb R $. This implies that if $ a < 0 $, we can choose $ x $ so that $ a < x < 0 $, and thus have $ g ( x ) = h ( x ) = 0 $, which gives $ a = b $. Similarly, if $ b < 0 $, choosing $ x $ with $ 0 < x < - b $ we get $ a = b $. Therefore, if $ a $ is not equal to $ b $ then $ a , b \ge 0 $.
By definition of $ g $ and $ h $ and the previous results, we have $ f ( x ) = a x $ for $ x \ge 0 $ and $ f ( x ) = b x $ for $ x < 0 $. If $ f $ is not linear, then $ a $ is not equal to $ b $, and thus we must have $ a , b \ge 0 $. We already know that whatever the value of $ a $ and $ b $, we'll get a solution in such a case. Thus there is no solution other than linear ones and those of this form, and we're done.
Best Answer
There is an easy to check polynomial solution. For $k=2^ho$ as above,
$$A(x):=\sum_{j=0}^h c^{h-j}x^{2^jo}=x^k+cx^{k/2}+c^2x^{k/4}+\dots+c^hx^{k/2^h}$$ Note that the only odd-degree term is the last one, $c^hx^o$.
$$\sim *\sim $$ If $\bf |c|<1$, this is also the only locally bounded solution on $\mathbb C$; more generally, for any $r\ge1$ and any bounded $\mathbb C$-valued function $f$ on the disk $B(0,r)$ of radius $r$, the linear functional equation $$A(x)=c\frac{A(\sqrt x)+A(-\sqrt x)}{2}+f(x)$$ has a unique bounded solution on $B(0,r)$ given inverting the operator $I-cH$ with $$(Hg)(x):=\frac{g(\sqrt x)+g(-\sqrt x)}{2}=\frac12\sum_{z^2=x}g(z), $$ the mean value of $g$ on the square roots of $x$. So the iterates of $H$ are given by $$H^n g(x)=\frac1{2^n}\sum_{z^{2^n}=x}g(z), $$ and give the mean value of $g$ on the $2^n$-th roots of $x$. Since $\|H\|\le 1$ one has the Neumann expansion $(I-cH)^{-1}=\sum_{n=0}^\infty H^n$ and finds $$A(x):=\sum_{n\in\mathbb N, z^{2^n}=x} (c/2)^n f(z).$$ If $f$ is a locally bounded function on $\mathbb C$, the latter is the only locally bounded solution on $\mathbb C$.