In Heaviside-Lorentz units, the equation for the electrostatic potential reads
\Delta \phi(\vec{x})=-\rho(\vec{x}).
In spherical coordinates, for a radially symmetric charge distribution, this equation simplifies to the ode,
\frac{1}{r} \mathrm{d}_r^2 (r \phi)=:\frac{1}{r} (r \phi)''=-\rho(r).
Here, r=|\vec{x}|. In your case we have
(r \phi)''=-k r^2.
This you can integrate up successively:
(r \phi)'=-\frac{k}{3} r^3+C_1, \quad r \phi=-\frac{k}{12} r^4 + C_1 r+C_2
or, finally,
\phi(r)=-\frac{k}{12} r^3 + C_1 + \frac{C_2}{r}.
The integration constants, C_1 and C_2, are determined by appropriate boundary conditions, where C_1 is physically irrelevant since the electric field is the observable quantity, which is given by the gradient of the potential. So I set C_1=0 Further, since the charge distribution is continuous at the origin, we must have C_2=0. So finally we get
\phi(r)=-\frac{k}{12} r^3, \quad \vec{E}=-\vec{\nabla} \phi = \frac{k}{4} r^2 \vec{e}_r = \frac{k}{4} r \vec{x}.
To check this result, we take the divergence, which for our case of \vec{E}=E_r(r) \vec{e}_r simplifies to
\vec{\nabla} \cdot \vec{E}=\frac{1}{r^2} (r^2 E_r)'=\frac{1}{r^2} \left (\frac{k r^4}{4} \right )'=k r=\rho
as it should be since Maxwell's equations for the electrostatic case read
\vec{\nabla} \times \vec{E}=0, \quad \vec{\nabla} \cdot \vec{E}=\rho.
I guess, your confusion comes from forgetting that the differential operators look different in curvilinear coordinates compared to the expressions in cartesian ones.
Of course, you can calculate the divergence also in Cartesian coordinates. Although this is a bit more inconvenient in your case of spherical symmetry, you must get the same result. Indeed using the Cartesian form
\vec{E}=\frac{k}{4} \sqrt{x^2+y^2+z^2} \vec{x}
you get again
\vec{\nabla} \cdot \vec{E}=\partial_x E_x+\partial_y E_y+\partial_z E_z=\frac{k}{4} \left (\frac{2 x^2+y^2+z^2}{r}+\frac{x^2+2y^2+z^2}{r}+\frac{x^2+y^2+2 z^2}{r} \right) =\frac{k}{4} \frac{4(x^2+y^2+z^2)}{r}=k r.
