In Heaviside-Lorentz units, the equation for the electrostatic potential reads
[tex]\Delta \phi(\vec{x})=-\rho(\vec{x}).[/tex]
In spherical coordinates, for a radially symmetric charge distribution, this equation simplifies to the ode,
[tex]\frac{1}{r} \mathrm{d}_r^2 (r \phi)=:\frac{1}{r} (r \phi)''=-\rho(r).[/tex]
Here, [itex]r=|\vec{x}|[/itex]. In your case we have
[tex](r \phi)''=-k r^2.[/tex]
This you can integrate up successively:
[tex](r \phi)'=-\frac{k}{3} r^3+C_1, \quad r \phi=-\frac{k}{12} r^4 + C_1 r+C_2[/tex]
or, finally,
[tex]\phi(r)=-\frac{k}{12} r^3 + C_1 + \frac{C_2}{r}.[/tex]
The integration constants, [itex]C_1[/itex] and [itex]C_2[/itex], are determined by appropriate boundary conditions, where [itex]C_1[/itex] is physically irrelevant since the electric field is the observable quantity, which is given by the gradient of the potential. So I set [itex]C_1=0[/itex] Further, since the charge distribution is continuous at the origin, we must have [itex]C_2=0[/itex]. So finally we get
[tex]\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}.[/tex]
To check this result, we take the divergence, which for our case of [itex]\vec{E}=E_r(r) \vec{e}_r[/itex] simplifies to
[tex]\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[/tex]
as it should be since Maxwell's equations for the electrostatic case read
[tex]\vec{\nabla} \times \vec{E}=0, \quad \vec{\nabla} \cdot \vec{E}=\rho.[/tex]
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
[tex]\vec{E}=\frac{k}{4} \sqrt{x^2+y^2+z^2} \vec{x}[/tex]
you get again
[tex]\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.[/tex]