In the context of this question I have considered doing some further calculations to check if I get somewhere. The first step was to solve the PDE
$$\frac{\partial \psi}{\partial z} = \frac{i}{2k} \frac{\partial^2 \psi}{\partial x^2}$$ which formally is a Schrödinger equation where the potential is $V(x) = 0$ and instead of time I have some spatial variable. This version comes from optics, but has a similar form.
At some point during my calculations I realized that I have to compute the solution of this PDE for an initial condition $\psi_0(x) = \delta(x)$ so I started doing that. My approach is based on using the Fourier Transform along the $x$-axis such that the PDE becomes
$$\frac{\partial \tilde{\psi}}{\partial z} = -\frac{i2\pi^2\xi^2}{k} \tilde{\psi}$$ which has the solution
$$\tilde{\psi} = \tilde{\psi_0} \exp \left(-\frac{i2\pi^2\xi^2}{k} z \right).$$
From this point the solution is retrieved using the inverse Fourier Transform which formally gives
$$\psi = \psi_0 * \mathscr{F}^{-1} \exp\left( -\frac{i2\pi^2\xi^2}{k} z\right)$$
where $*$ stands for convolution.
I have computed the inverse of the second term by hand and I have also checked it using Wolfram Alpha, and, as far as I can say, it should allow for the solution to be written as
$$\psi = \psi_0 * \left(\sqrt{\frac{k}{i2\pi z}}\exp (\frac{ikx^2}{2z})\right)$$
Now, if the initial conditions is $\psi_0(x) = \delta(x)$ as I have mentioned above, this means that the solution takes the form
$$\psi = \sqrt{\frac{k}{i2\pi z}}\exp \left(\frac{ikx^2}{2z}\right)$$
Now that I got the answer I looked at it and things did not look right. First of all, Why is the solution's amplitude decaying on the entire spatial domain proportional to $1/\sqrt{z}$?
Second, I did the computation for a gaussian initial condition with some width and if I take the width to be smaller and smaller, the solution spreads more rapidly, so I would assume that for a width that goes to 0 the solution should spread "instantly". But this implies that if I "measure" the position of a particle, got the delta Dirac function for the position and then compute how it evolves from that point onward, the position of my particle would be unknown and equally probable to be anywhere on the $x$-axis $\forall z>0$.
I am aware that my equation uses $z$ instead of $t$ and the constants are different, the reason being that my PDE comes from optics (paraxial approximation of the wave equation), but formally it is the same PDE with the same interpretation as for the quantum case, just scaled differently.
My guess here is that this method of solving the PDE fails when I do either the inverse Fourier Transform, either the convolution, because the second term is not square-integrable or something related to that. I am not sure however.
EDIT: In order to clarify the difference between my equation and the Schrödinger equation, since I mentioned that they are formally identical, I will write both of them as they are usually found. For Schrödinger with only one spatial coordinate and no potential term (I did simplify by $ih$):
$$\frac{\partial \psi}{\partial t} = i\frac{h}{2m} \frac{\partial^2 \psi}{\partial x^2}$$
My equation from optics:
$$\frac{\partial \psi}{\partial z} = \frac{i}{2k} \frac{\partial^2 \psi}{\partial x^2}$$
Now, ignoring the values of the constants, the two equations are identical in the sense that for the standard quantum picture $t$ is time and $x$ is space and the equation describes how the wavefunction evolves in time (along the $t$ axis), while for my case (optics) $z$ is the propagation axis and $x$ is the transverse axis and the equation describes how the optical wave evolves during propagation along the $z$ axis. So, if one solves one of the two forms, the solution should be, up to some scaling, the same for the other. This is the reason why I can consider the PDE as a IVP, not a BVP.