The complete solution for the time dependent equation with an infinite potential step is found by the method of images. Given any initial wavefunction
$$ \psi_0(x) $$
for x<0, you write down the antisymmetric extension of the wavefunction
$$ \psi_0(x) = \psi_0(x) - \psi_0(-x) $$
And you solve the free Schrodinger equation. So any solution of the free Schrodinger equation gives a solution for the infinite potential step. This is not completely trivial to make, because the solutions do not vanish in any region. But, for example, the spreading delta-function
$$ \psi(x,t) = {1\over \sqrt{2\pi it}} e^{-(x-x_0)^2\over it } $$
Turns into the spreading, reflecting, delta function
$$ \psi(x,t) = {1\over \sqrt{2\pi it}} e^{-(x-x_0)^2\over it } - {1\over \sqrt{2\pi it}} e^{-(x+x_0)^2\over it } $$
You can do the same thing with the spreading Gaussian wavepacket, just subtract the solution translated to +x from the solution translated to -x. In this case, normalizing the wavefunction is hard when the wavefunction start out close to the reflection wall.
Time independent infinite potential wall
The solution to the time independent problem of the infinite potential wall are all wavefunctions of the form
$$ \sin(kx) $$
for all k>0. Superposing these solutions gives all antisymmetric functions on the real line.
To find this solution, note that the time independent problem (eigenvalue problem) for the Schrodinger equation is solved by sinusoidal waves of the form $e^{ikx}$, and you need to superposes these so that they are zero at the origin, to obey the reflection condition. This requires that you add two k-waves up with opposite signs of k and opposite sign coefficients.
The opposite sign of k just means that the wave bounces off the wall (so that k changes sign), while the opposite sign of the coefficient means that the phase is opposite upon reflection, so that the wave at the wall cancels.
General solution
The time dependent problem for a time independent potnetial is just the sum of the solutions to the time independent problem with coefficients that vary in time sinusoidally.
If the eigenfunctions $\psi_n$ are known, and their energies $E_n$ are known, and the potential doesn't change in time, then the,
$$ \psi(t) = \sum_n C_n e^{-iE_n t} \psi_n(x) $$
is the general solution of the time dependent problem. This is so well known that generally people don't bother saying they solved the time-dependent problem once they have solved the eigenvalue problem.
The general solution of the time-dependent Schrodinger equation for time dependent potentials doesn't reduce to an eigenvalue problem, so it is a different sort of thing. this is generally what people understand when you say solving the time-dependent equation, and this reflects the other answers you are getting. I don't think this was the intent of your question, you just wanted to know how to solve the time dependent equation for a time independent potential, in particular, for an infinite reflecting potential wall. This is just the bouncing solution described above.
Note: I have tried to make my answer a little more general, with detail, so that it will be useful for more people.
The question is what boundary conditions do we apply to our wavefunction either side of a Dirac delta function?
In your example we have the potential $$V(x)=\begin{cases}\infty &\text{ if } x < 0\\ \alpha~\delta(x-a) &\text{ if } x \geq 0 \end{cases}$$ We are interested in the boundary conditions either side of $x=a$. What information do we have? Well, due to the probabilistic interpretation of the wavefunction we require continuity of the wavefunction. That is, our first boundary condition is
$$(1)~~~~~~~~~~~~\boxed{\psi_{-}(a) = \psi_{+}(a)}$$
where the $\pm$ subscripts represent the right and left sides of $x=a$ respectively. What other conditions can we set? Well, usually we would ask that the first derivative is also matched either side of $x=a$ (you should be asking yourself why do we do this?), but in this case this is not the right condition to impose. Let's see why.
Where does the boundary condition on $\frac{\partial \psi}{\partial x}$ come from?
Our wavefunction is a solution to the 1-D time-independent Schrödinger Equation:
$$H~\psi(x) = -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\psi(x)+V(x)\psi(x) = E~\psi(x)$$
Looking at this equation, we see that we can get a boundary condition on $\frac{\partial \psi}{\partial x}$ at any point $a$ by integrating it w.r.t $x$ over the region $[a-\epsilon,a+\epsilon]$, taking $\epsilon\rightarrow 0$:
$$\lim_{\epsilon\rightarrow 0}\left[-\int^{a+\epsilon}_{a-\epsilon}dx\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\psi(x)+\int^{a+\epsilon}_{a-\epsilon}dx~V(x)\psi(x)\right] = E~\lim_{\epsilon\rightarrow 0}~\int^{a+\epsilon}_{a-\epsilon}dx~\psi(x)\\
\implies -\frac{\hbar^2}{2m}\left[\frac{\partial \psi_{+}(a)}{\partial x}-\frac{\partial \psi_{-}(a)}{\partial x}\right] + \lim_{\epsilon\rightarrow 0}\int^{a+\epsilon}_{a-\epsilon}dx~V(x)\psi(x) = 0$$
where we have used the continuity of $\psi$ to evaluate the RHS as zero. Note that for any $V(x)$ the second term doesn't necessarily vanish. But, if $V(x)$ is continuous, this term will vanish for the same reason the RHS did, and in those cases we yield the usual boundary condition $$\boxed{\frac{\partial \psi_{+}(a)}{\partial x}=\frac{\partial \psi_{-}(a)}{\partial x}}~~~~~~~(\mbox{when $V(a)$ is continuous at $x=a$})$$
We note that the for example given here, $V(x)$ is definitely not continuous at $x=a$, where it diverges to infinity. So this second term doesn't vanish, in fact
$$\lim_{\epsilon\rightarrow 0}\int^{a+\epsilon}_{a-\epsilon}dx~V(x)\psi(x) = \alpha\lim_{\epsilon\rightarrow 0}\int^{a+\epsilon}_{a-\epsilon}dx~\delta(x-a)\psi(x) = \alpha~\lim_{\epsilon\rightarrow 0}~\psi(a) = \alpha~\psi(a)$$
so rearranging our results, the second boundary condition for the problem is
$$(2)~~~~~~~~~~~~\boxed{\frac{\partial \psi_{+}(a)}{\partial x}-\frac{\partial \psi_{-}(a)}{\partial x} = \frac{2m\alpha}{\hbar^2}\psi(a)}$$
Explicitly using the form of your wavefunctions (using $B=-A$ to eliminate $B$)
$$\psi_{+}(x) = Fe^{ikx},~~~~~~\psi_{-}(x) = A\left(e^{ikx}-e^{-ikx}\right)$$
this boundary condition becomes
$$\boxed{ikFe^{ika}-ik\left(Ae^{ika}+e^{-ika}\right) = \frac{2m\alpha}{\hbar^2}\psi(a)}$$
as required.
Best Answer
Yes, your model can be fully solved analytically, though actually finding the spectrum may involve the numerical solution of a trascendental equation. This is done explicitly in several introductory QM textbooks, but in essence all you need to do is solve locally for each region (which gives you sines and cosines, or exponentials) and then stitch the solutions by demanding continuity of the wavefunction and its derivative.
In essence, as $a$ becomes bigger than the ground state energy, you end up solving for two weakly coupled wells, with roughly independent ground states; the global well ground and first excited state then become even and odd linear combinations of the two individual ground states, with an energy splitting that depends on the coupling between the two wells.
Regarding your other points,
When the middle bit goes up to infinity.
When the middle bit goes up to infinity. The ground state is never zero there (though it does decay exponentially in $a$); the first excited state is identically zero in the middle.