Let wave function $\Psi$ be defined on domain $D \in \mathbb{R}^n$. The Neumann condition $\frac{\partial \Psi} {\partial {\bf n}} = 0$ on the boundary $\partial D$ has a simple interpretation in terms of the probability current of $\Psi$. For $\Delta \Psi = i \partial\Psi/\partial t$ (although it's usually taken as $i \partial\Psi/\partial t = - \Delta \Psi$), the probability current at an arbitrary point ${\bf x} \in \mathbb{R}^n$ is
$$
{\bf j}({\bf x}) = i [ \Psi^*({\bf x}){\bf \nabla}\Psi({\bf x}) - \Psi({\bf x}){\bf \nabla}\Psi^*({\bf x}) ]
$$
and the normal current on $\partial D$ reads
$$
{\bf n} \cdot {\bf j} = i\; [ \Psi^* \frac{\partial \Psi}{\partial {\bf n}} - \Psi \frac{\partial \Psi^*}{\partial {\bf n}} ]
$$
(has the wrong sign, I know, but I accounted for OP's form of the Sch.eq. as $\Delta \Psi = i \partial\Psi/\partial t$).
Setting $\frac{\partial \Psi} {\partial {\bf n}} = 0$ amounts to ${\bf n} \cdot {\bf j} = 0$ everywhere on $\partial D$, thus confining the corresponding system within $D$ without an infinite potential well, as under Dirichlet conditions ($\Psi = 0$ on $\partial D$). This is the case of perfect reflection on $\partial D$.
There is a mention of this in Section 5.2 of Visual Quantum Mechanics:
Selected Topics with Computer-Generated Animations of Quantum-Mechanical Phenomena, by Bernd Thaller (Springer, 2000); Google Books link.
As for applications, one answer to another post, Can we impose a boundary condition on the derivative of the wavefunction through the physical assumptions?, pointed to the use of Neumann conditions in R-matrix scattering theory.
$$
$$
Clarification following @arivero's observation on conditions necessary to trap the system within domain D:
We can say that the system described by $\Psi$ is trapped in domain D if the total probability to locate it in D, $P_D$, is conserved in time: $dP_D/dt = 0$. In this case, if the system is initially located within D, such that $\Psi({\bf x},t=0) = 0$ for all ${\bf x} \notin D$ and $P_D(t=0) = 1$, then it will remain in D at all $t > 0$, since $P_D(t) = P_D(t=0) = 1$. If initially $P_D(t=0) < 1$ (the system has a nonzero probability to be located outside D), then we still have $P_D(t) = P_D(t=0) < 1$.
Conservation of $P_D$ is equivalent to a condition of null total probability current through the boundary $\partial D$. Note that it is not necessary to require null probability current at every point of $\partial D$, but only null total probability current through $\partial D$.
The difference can be understood in terms of path amplitudes (path-integral representation). In the former case, the amplitude $\Psi({\bf x}_1 t_1, {\bf x}_2 t_2; {\bf x} t)$ that the system "goes" from point ${\bf x}_1 \in D$ at time $t_1$ to point ${\bf x}_2 \notin D$ at time $t_2 > t_1$ while passing through point ${\bf x} \in \partial D$ at some time $t$, $t_1 < t < t_2$, is nonzero $\forall {\bf x} \in \partial D$: $\Psi({\bf x_1} t_1, {\bf x_2} t_2; {\bf x} t)≠0$. If however we demand null probability current at every point of $\partial D$, then $\Psi({\bf x}_1 t_1, {\bf x}_2 t_2; {\bf x} t)=0$, $\forall {\bf x} \in \partial D$.
In other words, null total probability current on $\partial D$ enforces weak trapping in the sense that overall $P_D(t) = $ const. and "cross-over events" across the boundary balance out. Null local probability current at every point of $\partial D$, ${\bf n} \cdot {\bf j} = 0$, corresponds to strong trapping in the sense that the system "does not cross" at any point ${\bf x} \in \partial D$. Imposing the strong trapping condition is equivalent to requiring that the weak trapping condition be satisfied by any wave function $\Psi$, as opposed to one selected $\Psi$. In this case the system is essentially confined within D at all times. Incidentally, the strong trapping condition follows from the requirement that the restriction of the system Hamiltonian on domain D remain self-adjoint.
Derivation of probability current conditions:
The free Schroedinger equation for $\Psi$, $i\partial\Psi/\partial t = \Delta\Psi$ as above (OP's choice of sign), implies local conservation of the probability density $\rho({\bf x,t}) = \Psi({\bf x},t)\Psi^*({\bf x},t)$:
$$
\frac{\partial \rho({\bf x},t)}{\partial t} + {\bf \nabla}\cdot {\bf j}({\bf x},t) = 0
$$
Integrating this over domain D yields
$$
\int_D dV\;\frac{\partial {\rho({\bf x},t)}}{\partial t} + \oint_{\partial D}dS\;{{\bf n}\cdot{\bf j}} = \frac{d}{dt}\int_DdV\;{\rho({\bf x},t)} + \oint_{\partial D}dS\;{{\bf n}\cdot{\bf j}} = 0
$$
which after denoting $P_D = \int_DdV\;{\rho({\bf x},t)}$ becomes
$$
\frac{dP_D}{dt} + \oint_{\partial D}dS\;{{\bf n}\cdot{\bf j}} = 0
$$
Imposing $dP_D/dt = 0$ necessarily means $\oint_{\partial D}dS\;{{\bf n}\cdot{\bf j}} = 0$. Note that $\oint_{\partial D} dS\;{{\bf n}\cdot{\bf j}} = 0$ does not require ${{\bf n}\cdot{\bf j}} = 0$ at every point on $\partial D$, whereas ${{\bf n}\cdot{\bf j}} = 0$ does imply $\oint_{\partial D} dS\;{{\bf n}\cdot{\bf j}} = 0$ and $dP_D/dt = 0$.
Self-adjoint restriction of the free Hamiltonian on domain D:
A self-adjoint restriction of $H\Psi = \Delta \Psi$ on D requires that $
\int_D dV\;\Phi^* (\Delta\Psi) = \int_D dV\;(\Delta\Phi^*) \Psi$ or $
\int_D dV\;[\Phi^* (\Delta\Psi) - (\Delta\Phi)^* \Psi] = 0$.
Use $\Phi^* (\Delta\Psi) = {\bf \nabla}\cdot(\Phi^* {\bf \nabla}\Psi) - {\bf \nabla}\Phi^* \cdot {\bf \nabla}\Psi$ and Green's theorem to obtain
$$
\int_D dV\;[\Phi^* (\Delta\Psi) - (\Delta\Phi)^* \Psi] = \oint_{\partial D} dS\;[\Phi^*\frac{d\Psi}{d{\bf n}} - \Psi\frac{d\Phi^*}{d{\bf n}}] = 0
$$
If the last condition above is to be satisfied by arbitrary $\Phi$, $\Psi$, it must hold locally:
$$
\Phi^* \frac{d \Psi}{d{\bf n}} - \Psi \frac{d \Phi^*}{d{\bf n}} = 0
$$
This means $\frac{1}{\Psi}\frac{d\Psi}{d{\bf n}} = \frac{1}{\Phi^*}\frac{d\Phi^*}{d{\bf n}} = a({\bf x}), \forall {\bf x} \in \partial D$. The case $\Phi = \Psi$ shows that $a({\bf x}) = a^*({\bf x})$. Therefore the restriction of $H$ on $D$ is self-adjoint if and only if wave functions in its support satisfy a boundary condition
$$
\frac{d\Psi}{d{\bf n}} = a({\bf x})\Psi
$$
for some given real-valued function $a({\bf x})$. In particular, this means every wave function $\Psi$ also satisfies the strong trapping condition ${\bf n} \cdot {\bf j} = 0$.
Finally note that the strong trapping condition means $\Psi^*\frac{d\Psi}{d{\bf n}} - \Psi\frac{d\Psi^*}{d{\bf n}} = 0$, but does not imply that $\Phi^*\frac{d\Psi}{d{\bf n}} - \Psi\frac{d\Phi^*}{d{\bf n}} = 0$ for arbitrary $\Psi$, $\Phi$. If we consider the latter expression as the matrix element of a "local current operator" ${\bf n} \cdot {\bf \hat j}$, then the strong trapping condition requires that diagonal elements of ${\bf n} \cdot {\bf \hat j}$ are 0, whereas self-adjointness requires that the entire space of wave functions is in the kernel of each ${\bf n} \cdot {\bf \hat j}$.
Best Answer
I'm not finding solving this enormously difficult but it is both tedious and very lengthy. For that reason I'll take the unusual approach to publish the solution in 3 parts:
Part 1: up to solving the homogeneous PDE.
Part 2: solving the non-homogeneous PDE.
Part 3: testing the solution with a simple example.
Problem statement: $$u_t=ku_{xx}$$ With boundary and in initial conditions: $$u(x,0)=f(x)$$ $$u(0,t)=b_1(t),u(L,t)=b_2(t)$$ So we're looking to solve the heat equation in one dimension, without heat sinks or sources but with time-dependent boundary conditions.
Part 1:
We assume the solution to be of the partitioned form: $$u(x,t)=s(x,t)+v(x,t)$$ So that: $$s_t+v_t=ks_{xx}+kv_{xx}...(i)$$ Translating the boundary/initial conditions: $$s(x,0)+v(x,0)=f(x)$$ $$s(0,t)+v(0,t)=b_1(t)$$ $$s(L,t)+v(L,t)=b_2(t)$$ We assume $s(x,t)$ to be linear in $x$ but time-dependent, so the general form of $s(x,t)$ is: $$s(x,t)=a(t)x+b(t)$$ So: $$s_{xx}=0$$ So with $s(0,t)$ and $s(L,t)$ we get: $$b(t)=b_1(t)$$ $$a(t)L+b(t)=b_2(t)$$ $$a(t)=\frac{b_2(t)-b_1(t)}{L}$$ So: $$s(x,t)=\frac{b_2(t)-b_1(t)}{L}x+b_1(t)$$ With $(i)$ and $s_{xx}=0$, we get: $$v_t=kv_{xx}-s_t$$ We solve the homogeneous equation first: $$v_t=kv_{xx}$$ With boundary conditions: $$v(0,t)=0,v(L,t)=0...(ii)$$ Using the Ansatz: $$v(x,t)=X(x)T(t)$$ Obtaining separation of variables: $$\frac1k \frac{T'}{T}=\frac{X''}{X}=-m^2$$ ODE for $X(x)$: $$X''+m^2X=0$$ $$X=c_1\cos mx+c_2\sin mx$$ With boundary conditions $(ii)$: $$c_1=0$$ $$m=\frac{n\pi}{L}, n=1,2,3....$$ $$X_n(x)=\sin\Big(\frac{n\pi x}{L}\Big)$$
Part 2:
Now for the non-homogeneous PDE: $$v_t=kv_{xx}-s_t$$ With boundary/initial conditions: $$v(0,t)=0,v(L,t)=0$$ $$v(x,0)=f(x)-s(x,0)$$ Assume the solution to be of the form: $$v(x,t)=\displaystyle \sum_{n=1}^{+\infty}T_n(t)X_n(x)$$ $$s_t=\frac{b'_2(t)-b'_1(t)}{L}x+b'_1(t)=-\displaystyle \sum_{n=1}^{+\infty}Q_n(t)X_n(x)$$ $$Q_n(t)=-\frac{\int_0^Ls_tX_n(x)dx}{\int_0^LX_n(x)X_n(x)dx}$$ Our solution must of course obey the original equation, so: $$\frac{\partial}{\partial t}\displaystyle \sum_{n=1}^{+\infty}T_n(t)X_n(x)=k\frac{\partial^2}{\partial t^2}\Bigg[\displaystyle \sum_{n=1}^{+\infty}T_n(t)X_n(x)\Bigg]+\displaystyle \sum_{n=1}^{+\infty}Q_n(t)X_n(x)$$ We know that: $$X''_n(x)=-m^2X_n(x)$$ So: $$\displaystyle \sum_{n=1}^{+\infty}T'_n(t)X_n(x)=\displaystyle \sum_{n=1}^{+\infty}-km^2T_n(t)X_n(x)+\displaystyle \sum_{n=1}^{+\infty}Q_n(t)X_n(x)$$ $$\displaystyle \sum_{n=1}^{+\infty}[T'_n(t)+km^2T_n(t)]X_n(x)]=\displaystyle \sum_{n=1}^{+\infty}Q_n(t)X_n(x)$$ So that: $$T'_n(t)+km^2T_n(t)=Q_n(t)$$ Solve with an integration factor to: $$T_n(t)=e^{-km^2t}\int_0^te^{km^2t}Q_n(t)dt+C_ne^{-km^2t}$$ Using the initial condition: $$v(x,0)=f(x)-s(x,0)=\displaystyle \sum_{n=1}^{+\infty}T_n(0)X_n(x)=\displaystyle \sum_{n=1}^{+\infty}C_nX_n(x)$$ So: $$C_n=\frac{\int_0^L[f(x)-s(x,0)]X_n(x)dx}{\int_0^LX_n(x)X_n(x)dx}$$
Part 3:
In Part 3 I'll test drive the solution for the simple case:
$$u(x,0)=u_0$$ $$u(0,t)=u_0-at$$ $$u(L,t)=u_0+bt$$
Bits and pieces we need: $$s(x,t)=\frac{b_2(t)-b_1(t)}{L}x+b_1(t)=u_0+\frac{a+b}{L}xt-at$$ $$s_t=\frac{b'_2(t)-b'_1(t)}{L}x+b'_1(t)=\frac{a+b}{L}x-a$$ $$s(x,0)=u_0$$ $$Q_n(t)=-\frac{\int_0^Ls_tX_n(x)dx}{\int_0^LX_n(x)X_n(x)dx}$$ Because the boundary conditions are only linear in $t$ the $Q_n$ turn out time independent and compute to: $$Q_n=2\Big[\frac{b}{n\pi}\cos n\pi+\frac{a}{n\pi}\Big]$$ Because: $$s(x,0)=u_0=f(x)\implies C_n=0$$ So slightly reworked we get: $$T_n(t)=\frac{2L^2}{kn^3\pi^3}\Big[(-1)^nb+a\Big]\Big[1-e^{-k\Big(\frac{n\pi}{L}\Big)^2t}\Big]$$
Put it altogether with: $$u(x,t)=s(x,t)+\displaystyle \sum_{n=1}^{+\infty}T_n(t)X_n(x)$$ Where $X_n(x)=\sin\Big(\frac{n\pi x}{L}\Big)$
Let's evaluate for $u_0=100$, $L=1$, $k=1$, $a=1$ and $b=2$.
Firstly let's look at $s(x,t)$. The significance of $s(x,t)$ is as follows. If we allow temperature evolution of the boundaries for some time $t$, then stop heating/cooling and maintain the boundary temperatures until thermal equilibrium is achieved, then $s(x,t)$ describes the temperature profile, e.g.:
Here's $u(x,t)$ with the first three terms computed:
And here's $v(x,t)$ with the first three terms computed: