You have received other answers, so I would like to focus on the overall problem of necessary regularity of the solutions of Schroedinger equation. More precisely: why one should require the regularity conditions on $\psi$ you have read in the other answers?
The point can be traced back to one of the most important axioms of QM: observables are self-adjoint operators. The reason for this requirement is that self-adjoint observables admits a spectral decomposition in terms of orthogonal projectors labeled by (Borelian) subsets of $R$ interpreted as the set of results of the measurement of the observable.
(It is possible to make weaker the requirement dealing with decompositions of bounded positive operators, but I only stick to the elementary case here.)
In our case the relevant observable is the Hamiltonian one. But for the sake of simplicity I intend to start focusing on the momentum observable along the $k$th axis: $M_k$. One usually assumes that:
$$M_k := -i\hbar \frac{\partial}{\partial x_k}\:,\qquad (0)$$
where, for instance the domain ${\cal D}(M_k)$ is $C_0^\infty(R^3)$ of $\cal S(R^3)$ (the Schwartz'space), what follows is independent on this choice.
It is true that, if $\psi,\phi \in {\cal D}(M_k)$ then:
$$\langle \psi| M_k \phi \rangle = \langle M_k\psi | \phi \rangle\:.$$
Indeed, that identity only says that $M_k$ is symmetric (a densely defined operator is symmetric if, on its domain, coincides to the adjoint operator). However, it does not say that $M_k$ is self-adjoint. The self-adjoint condition (that implying the existence of the spectral decomposition) is instead:
$$M_k^\dagger = M_k\:. \qquad (1)$$
Above $M^\dagger$ is defined as, follows. First one defines its domain:
$${\cal D}(M^\dagger_k) := \left\{\psi \left.\in L^2(R^3) \:\right|\: \exists \psi' \in L^2(R^3)\qquad \mbox{with}\quad\langle \psi' | M_k \phi\rangle = \langle \psi |\phi\rangle\quad \forall \phi \in {\cal D}(M_k)\right\}$$
Since ${\cal D}(M_k)$ is dense, $\psi'$ is uniquely determined by $\psi$ and thus the map:
$${\cal D}(M_k^\dagger) \ni \psi \mapsto \psi' =: M_k^\dagger \psi$$
is well-defined. One easily sees that ${\cal D}(M_k^\dagger)$ is a subspace of $L^2(R^3)$ with ${\cal D}(M_k^\dagger) \supset {\cal D}(M_k)$ and that
$M_k^\dagger$ is linear.
In this case ${\cal D}(M_k^\dagger)$ turns out to be considerably greater than ${\cal D}(M_k)$, so that (1) fails and $M_k$ is not self-adjoint with the given (standard) definition. What it is true is that $M_k^\dagger$, defined as above, is self-adjoint and that it is the only self-adjoint extension of $M_k$.
Mathematically one says that $A$ is essentially self-adjoint when it is symmetric its adjoint operator $A^\dagger$ is self-adjoint $A^\dagger = (A^\dagger)^\dagger$. Therefore $M_k$ is essentially self-adjoint.
This discussion leads to conclude that the true definition of the momentum operator is not (0) but is:
$$P_k = M^\dagger_k\:.$$
It is however important to stress that the naive, technically wrong, definition (0) uniquely defines $P_k$, since it is the only self-adjoint extension of $M_k$. This latter is very simple to handle, because is a differential operator. Conversely $P_k$ has a domain much more difficult to characterize (without using the Fourier transform). The domain of $P_k$ is made of the functions $\psi$ in $L^2(R^3)$ that admit weak $k$-derivative which, in turn, is a function in $L^2(R^3)$.
One says that $\psi \in L^2(R^3)$ admits a weak $k$-derivative $\phi_k : R^3 \to C$, if there is a function - the mentioned $\phi_k$ - such that
$$\int_{R^3} \frac{\partial f}{\partial x_k} \psi \:d^3x = - \int_{R^3} f \phi_k \:d^3x \quad \forall f \in C_0^\infty(R^3)\:.$$
You see that if $\psi$ admits the standard $k$-derivative it coincides with the weak one (that therefore exists in this case). However, there are many functions admitting weak derivatives that are not differentiable anywhere!
Let us come to the problem of the Hamiltonian operator. The Hamiltonian operator, in the mathematically naive version for the non relativistic theory, always includes an added part proportional to the Laplacian operator $\Delta$. As a matter of fact, for $a := -\hbar^2/(2m)$ and for some function $V: R^3 \to R$ the naive Hamiltonian operator is:
$$A := a \Delta + V\:,$$
with domain ${\cal D}(A)$ made of sufficiently differentiable functions.
Again one would be sure that $A$ is self-adjoint to exploit all spectral technology, but, as before $A$ is not. At most, with a careful choice of the domain ${\cal D}(A)$, the operator $A$ turns out to be essentially self adjoint. Namely,
$A^\dagger$ is self-adjoint and the true Hamiltonian observable can safely be defined as: $$H := A^\dagger\:.$$
As before, the correct domain (and one could have many choices!) involves weak derivatives: $\Delta$ has to be interpreted using (second) weak derivatives instead of standard derivatives. So the class of functions one should consider in solving problems like that of finding the eigenvalues of $H$ (the energies of stationary states) or other problems, like determining the scattering states, are a wide class of generally non differentiable functions.
This is not the whole story, because, differently from the case of the momentum operator, the presence of $\Delta$ in $A$ makes easier the problem in view of known results on elliptic regularity. The basic results (due to Weyl, Friedrichs and Sobolev) under suitable hypotheses establish that if a function (actually a distribution) in $R^n$ verifies an equation like $$\Delta f = g\:,$$
where $\Delta$ is interpreted in weak sense, then the degree of weak differentiable regularity of $f$ is that of $g$ plus $2$. Moreover, if $f$
(assumed to be locally $L^2$) has a certain degree $k$ of weak regularity, it also has another degree $k' = k-p$ of standard regularity, where $p>0$ is a number depending on $n$.
(To write a rigorous statement I should introduce several mathematical notions and I will not do for the sake of simplicity, since I just wish to give an idea about the basic argument).
Taking this result into account, it turns out that, for instance, if $\psi \in {\cal D}(H)= {\cal D}(A^\dagger)$ is an eigenvector of $H$, so that
$$-a\Delta_w \psi = (E-V) \psi \quad \mbox{where $\Delta_w$ is the weak Laplacian}\:,$$
then $\psi \in C^\infty$ where $V$ is such.
These procedures and results lead to a precise theorem concerning the mathematical requirements
for a function $\psi$ that stays in the domain of $H$ and, if it the case, solves the proper or generalized eigenvalue equation. The theorem takes into account the fact that the true self-adjoint Hamiltonian operator is not the differential operator $A$, but is its unique self-adjoint extension $A^\dagger$.
The theorem considers an operator of the form:
$$A = a \Delta + V$$
where $V: R^3 \to R$ has the form for $N$ real constants $g_k$ and corresponding isolated points ${\bf x}_k$:
$$V({\bf x}) = \sum_{j=1}^N \frac{g_k}{|{\bf x}-{\bf x_j}|} + V_0({\bf r})\:,$$
$V_0$ is bounded below, diverges at most polynomially for $|{\bf x}|\to +\infty$ and it is a continuous function except for a finite number of 2-surfaces $\Sigma_i$ where the discontinuities are finite. With these hypotheses it is possible to establish that $A$ is essentially self adjoint on ${\cal D}(A)= C_0^\infty(R^3)$ or ${\cal D}(A)={\cal S}(R^3)$ with the same unique self-adjoint extension $H= A^\dagger$ in both cases. The domain of $H$ is much larger than these spaces and includes functions that so not admit proper second derivatives in the whole $R^3$.
It turns out that the functions (distributions in the generalize case) $\psi : R^3 \to C$ that can be used to solve the proper or generalized eigenvalue problem for the true Hamiltonian $H := A^\dagger$ must verify the following requirements (in addition to the eigenvalue problem):
1) away from the singularities of $V$, $\psi$ is $C^2$ and solves the eigenvalue equation interpreting $\Delta$as a proper differential operator;
2) crossing a singular surfaces $\Sigma_i$, if ${\bf y}\in \Sigma_i$, the function $\psi$ satisfies
$$\lim_{{\bf x} \to {\bf y}^+}\psi({\bf x})= \lim_{{\bf x} \to {\bf y}^-}\psi({\bf x})$$ where the two limits are computed from the two semi-spaces separated by $\Sigma$ around ${\bf y}$ and, similarly:
$$\lim_{{\bf x} \to {\bf y}^+}{\bf n}\cdot \nabla\psi({\bf x})= \lim_{{\bf x} \to {\bf y}^-}{\bf n}\cdot \nabla \psi({\bf x})$$
where ${\bf n}$ is the unit vector normal to $\Sigma_i$ at ${\bf y}$;
3) If ${\bf x}_k$ is an isolated singular point for $V$, the limit of $\psi$ for ${\bf x} \to {\bf x}_k$ exists and is finite.
Dealing with genuine $1D$ systems one finds similar requirements on the allowed wavefunctions.
Taking the above results into account, you can understand why, for instance a wavefunction in $R^3$ must not diverge at the origin: it is nothing but the requirement (3) above or even a consequence of requirement (1) when $V$ is regular at ${\bf r}=0$. For this reason $\psi(r) \sim r^{-1}\cos (kr)$ for $r\to 0$ cannot be accepted even if the corresponding $1D$ wavefunction $r\psi(r)$ is, in principle, allowed. The comparison from the $1D$ and the $3D$ case, based on the replacement $\psi(r) \to r\psi(r)$ is only formal and can only be used outside the singularities, while what is permitted or forbidden crossing a singularity must be discussed separately, remembering the true nature of the problem: $3D$ or $1D$.
Notice also that, for the free particle, the accepted solution for $\ell=0$: $$\psi({\bf x})= A\frac{\sin (kr)}{r}$$
in accordance with (1) is $C^2$ (more strongly it is real analytic) and not only bounded in a neighborhood of ${\bf x}=0$, including that point. In fact there is no singularity in ${\bf x}=0$ for the free particle since $V\equiv 0$ in that case, and the apparent singularity is only due to the use of polar coordinates.
