Something wrong with solving the Schrödinger equation 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.
 A: Your calculations are all correct (from what I can tell) and completely applicable. Moreover, your conclusion is also correct:

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$.

This is true, and it is simply an expression of the Heisenberg uncertainty principle. If you measure the position with an apparatus with uncertainty $\Delta x$, you will collapse the state to a localized wavefunction of that size, and then the uncertainty principle requires the momentum uncertainty $\Delta p \geq \hbar/\Delta x$ to grow as $\Delta x$ gets smaller. Thus, if $\Delta x$ becomes small, then the particle has a large velocity spread, and it will diffuse very quickly after measurement to a very large spatial spread.
Moreover, things are easier if you don't involve matter waves into any of this. If you come from optics, the uncertainty principle is entirely uncontroversial (you just substitute the wavenumber $k$ instead of the momentum, and forget about $\hbar$), and in this incarnation it simply describes single-slit diffraction: the narrower the slit, the faster the diffraction.

Edit in response to your comment,

Just to make sure, one thing to me is not as obvious. How does the decrease of the amplitude proportional to $1/\sqrt{z}$ fit in this picture? For the gaussian case I calculate the solution, apply a limit and, if I understand it correctly, I should get an amplitude that goes to 0 because it spreads on the entire domain. Here on the other hand it is already constant on the entire domain, but it starts from an infinite value and decays.

The state is not normalizable to begin with, so worrying about the norm is pretty pointless. That said, the local decrease of the amplitude does carry significance, not through the global norm but rather through the local conservation of probability (i.e. the continuity equation for $|\psi|^2$, which holds regardless of global normalizability as a consequence of the Schrödinger equation). The $1/\sqrt{z}$ dependence of $\psi$ tells you that $\nabla\cdot\vec j$ is nonzero, i.e. more particles leave any stretch of line than enter it from the other side.
This comes from the phase profile of the wavefunction (depicted in the figure here), which has the higher-momentum components increasingly far from the center, and diffusing outwards. The probability at each point is decreasing because that population is moving outwards, with not-quite-so-much population diffusing into that line element from the center. For a gaussian wavepacket, even in the flat center, this is fine as there is ultimately an edge that's expanding outwards.
For the delta spike, there is no edge, but that doesn't mean you can't have a conveyor belt transferring population out to infinity without end, not unlike the guests in Hilbert's hotel.
A: Your PDE:
$$\frac{\partial \psi}{\partial z} = \frac{i}{2k} \frac{\partial^2 \psi}{\partial x^2}$$
is solvable by separation of variables. Assume ('Ansatz'):
$$\psi(x,z)=X(x)Z(z)$$
Insert into the PDE:
$$XZ'=\frac{i}{2k}ZX''$$
Divide through by $XZ$:
$$\frac{Z'}{Z}=\frac{i}{2k}\frac{X''}{X}$$
$$\frac{2k}{i}\frac{Z'}{Z}=\frac{X''}{X}=-\lambda^2$$
where $\lambda^2$ is a separation constant. This gives us:
$$X''(x)+\lambda^2 X(x)=0\tag{1}$$
$$\frac{Z'(z)}{Z(z)}=-\lambda^2 \frac{i}{2k}\tag{2}$$
From $(1)$:
$$X(x)=A\sin\lambda x+B\cos\lambda x$$
And from $(2)$:
$$Z(z)=C\exp\Big(-\lambda^2 \frac{i}{2k}z\Big)$$
Now, to determine the integration constants $A$, $B$ and $C$, as well as the separation constant, we need boundary conditions. $2$ for $X(x)$ and $1$ for $Z(z)$ but these appear not to have been provided.
