3d solutions for 1d Schrödinger equation?

The general approach is that for Schrödinger equations where the potential is separable (in the sense that $V(x,y,z) = V_1(x) + V_2(y) + V_3(z)$) then there exists a basis of hamiltonian eigenfunctions which are separable (in the sense that $\psi(x,y,z) = \phi(x)\chi(y)\xi(z)$). However, in general, there are also non-separable eigenfunctions of the hamiltonian.

As regards the time-dependent Schrödinger equation, the details depend not only on the potential, but also on the initial condition. There are plenty of separable solutions, and if the initial condition is separable then the solution will remain separable. Conversely, if you start with a non-separable initial condition then the solution will remain non-separable.

The separability of the time-independent equation is handled in detail in every textbook so instead I'll show how this works for the time-dependent version. Suppose that we start with the Schrödinger equation in the form $$ i\hbar \frac{\partial}{\partial t}\psi(x,y,z,t) = \left[ -\frac{\hbar^2}{2m}\nabla^2 + V_1(x) + V_2(y) + V_3(z) \right]\psi(x,y,z,t) . \tag 1 $$ If you want a general solution to this equation, you need to specify an initial condition. In the absence of that, let's explore some particular solutions, and particularly, let's explore separable ones, i.e., solutions of the form $$ \psi(x,y,z,t) = \phi(x,t)\chi(y,t)\xi(z,t). \tag 2 $$ If you plug this into $(1)$, it's easy to see that a sufficient condition for $(1)$ to hold is if each of the individual 1D Schrödinger equations hold: \begin{align} i\hbar \frac{\partial}{\partial t}\phi(x,t) & = \left[ -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V_1(x)\right]\phi(x,t) \\ i\hbar \frac{\partial}{\partial t}\chi(y,t) & = \left[ -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial y^2} + V_2(y)\right]\chi(y,t) \tag 3 \\ i\hbar \frac{\partial}{\partial t}\xi(z,t) & = \left[ -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial z^2} + V_3(z)\right]\xi(z,t) . \end{align} (This also turns out to be a necessary condition. The full equation $(1)$, when divided by $\psi(x,y,z,t)$, comes down to a sum of three terms, each of which depends exclusively on $x$, $y$ and $z$, respectively, at fixed $t$. This is only possible if all three of the terms are uniformly zero.)

How does this relate to your question? In your example, $V_2(y)=0=V_3(z)$, so you can find a basis of TDSE solutions of the form $$ \chi_k(y,t)=e^{i(ky-\omega_k t)}, \quad \xi_k(z,t)=e^{i(kz-\omega_k t)}, $$ with $\omega_k = \frac{\hbar}{2m} k^2$. The specific example you've found uses the special case of $\chi_k(y,t)$ and $\xi_k(z,t)$ with $k=0$. This acts to mask what is really happening: your solution looks like a 1D problem, because it is actually three 1D solutions in tensor product with each other, with two of those being trivial.

So, with that as background, to address your question:

Does that mean that even for 1d potentials, one can have solutions which are not 1d?

yes, absolutely. Any solution of the $y$ and $z$ Schrödinger equations will work here.

Now, there is still a sense in which those solutions are "effectively 1D", though, in the sense that none of the separate 1D Schrödinger equations talk to each other, and the wavefunction remains separable. And this raises the question: are there any solutions which are not separable?

The answer there, again, is: yes, absolutely. Because of the linearity of the Schrödinger equation, given any two separable TDSE solutions $\psi_1(x,y,z,t) = \phi_1(x,t)\chi_1(y,t)\xi_1(z,t)$ and $\psi_2(x,y,z,t) = \phi_2(x,t)\chi_2(y,t)\xi_2(z,t)$, their linear combination $$ \psi(x,y,z,t) = \psi_1(x,y,z,t) + \psi_2(x,y,z,t) $$ is also a TDSE solution. And, as it turns out, if the individual components $\psi_1(x,y,z,t)$ and $\psi_2(x,y,z,t)$ are different enough (say, as one possible sufficient condition, $\chi_1(y,t)$ and $\chi_2(y,t)$ are orthogonal) then one can prove that the linear combination $\psi(x,y,z,t)$ cannot be written out as a product of individual 1D solutions.

The answer is no. You are mixing up domains of your functions, which is why you are getting such a result. There is a very large difference between a three-dimensional potential that depends only on $x$ and a proper one-dimensional potential. Recall that a function is defined by stating domains and then a rule. The rule may be the same, but the domains differ. For example, if we have $f(x)=x$ and $g(x,y,z)=x$, these are very different functions; one is a function from $\mathbb{R}$ to $\mathbb{R}$, while the other is a function from $\mathbb{R}^3$ to $\mathbb{R}$ ($f:\mathbb{R}\rightarrow\mathbb{R}$, $g:\mathbb{R}^3\rightarrow\mathbb{R}$).

Let's say we have two potentials, $V_1$ and $V_2$. These give the same results for all $x$, but $V_1$ is a function of one variable and $V_2$ a function of three. For the one-dimensional potential, we have $$ \begin{align} V_1&: \mathbb{R}\rightarrow\mathbb{R}\\ i\hbar\partial_t\langle x|\phi(t)\rangle&=-\frac{\hbar^2}{2m}\partial^2_x\langle x|\phi(t)\rangle+V_1(x)\langle x|\phi(t)\rangle \end{align} $$ Note that this is in one-dimension: the Laplacian becomes $\partial^2_x$ and $|\phi(t)\rangle$ is a function of $x$ alone. Now, let's move to three dimensions: $$ \begin{align} V_2&: \mathbb{R}^3\rightarrow\mathbb{R}\\ i\hbar\partial_t\langle x,y,z|\psi(t)\rangle&=-\frac{\hbar^2}{2m}\nabla^2\langle x,y,z|\psi(t)\rangle+V_2(x)\langle x,y,z|\psi(t)\rangle \end{align} $$ Note that $V(x)$ remains exactly the same in form, but it is now in the context of three-dimensional space; we've expanded the domain. Thus, now the solutions are of the form $\psi(x,y,z,t)$, but the solutions in the first case are in the form $\phi(x,t)$. Again, these are not the same, and $V_1$ and $V_2$ are not the same; while they give the same answer for any $x$, their domains are different and thus they are different functions. So no, you cannot have three-dimensional solutions to a one-dimensional problem.

To show this, we just need to show that $\partial\psi/\partial y$ and $\partial\psi/\partial z$ are zero.

No, this is wrong. These don't need to be zero.

Instead, you can solve the Schrödinger equation $$i\hbar\frac{\partial\psi}{\partial t}(x,y,z,t)= -\frac{\hbar^2}{2m}\nabla^2\psi(x,y,z,t)+\mathcal{V}(x)\psi(x,y,z,t).$$ by separation of variables with the approach $$\psi(x,y,z,t)=A(x)B(y)C(z)D(t).$$ where $A$, $B$, $C$ and $D$ are unknown functions of only one variable.
Then you easily find the solutions for the $y$-, $z$- and $t$-dependent parts $$\begin{align} B(y)&=B_0e^{ik_y y} \\ C(z)&=C_0e^{ik_z z} \\ D(t)&=D_0e^{-i\omega t} \end{align}$$ where $k_x$, $k_y$ and $\omega$ are arbitrary real constants.
And you are left with an ordinary differential equation for the $x$-dependent part: $$\left(\hbar\omega-\frac{\hbar^2(k_y^2+k_z^2)}{2m}\right)A(x)= -\frac{\hbar^2}{2m}\frac{d^2A(x)}{dx^2}+\mathcal{V}(x)A(x).$$