The solver uses a collocation method, that is, at each interval of the mesh it discretizes the differential equation by approximating the solution by a polynomial of a fixed degree and adding the necessary number of equations to equal the number of free parameters by demanding exact equality of the ODE at some interior points of the interval. For the given polynomial ODE this results in algebraic equations for the coefficients of the polynomial. The solver then proceeds to find a numerical solution of the BVP by solving this large non-linear system. After that, the error between the mesh points is estimated and the mesh refined towards a more homogeneous error estimate.
The point now is that the discretized differential equation will have more than one solution. Typically, the other solutions are not smooth, have large jumps between them, so that they are easily detectable by the solver. Or more often, the initial approximation is much closer to the smooth solution than to any of the extraneous solutions.
However, for stiff systems that might no longer be the case. "stiff" means that you have fast and slow components. The fast components require a very small step size to be modeled correctly in the discretization, while the many steps that that requires might lead to an undue accumulation of floating point errors in the other components. So it requires extra effort to find a workable balance.
In the collocation model, the large derivatives of the fast components lead to the situation where even the sampled exact solution is visually not smooth. That is, even some of the normal steps of the exact solution look like big jumps. So the true numerical solution is no longer easily distinguishable from the extraneous jumpy solutions of the algebraic system.
In the given system, the exact solution for large $K$ jumps at $x=1$ from essentially $c_2(x)\approx 0$ to $c_2(1)$ in the last part of the interval of length about $\Delta t=1/\sqrt{K}$. The square root is arbitrary for $1\gg\Delta t\gg1/K$, but captures the essence of the situation.
Close to $x=1$ (there can be no boundary layers at other points of the interval) the slow component $c_1$ can be treated as a constant, so that $c_2$ behaves qualitatively like an exponential function
$$
c_2(x)=c_2(1)\,e^{K\,c_1(1)\,(x-1)}=c_2(1)\,e^{-K\,c_1(1)\,(1-x)}
$$
This situation is captured by $e^{-Kt}$ close to $t=0$, i.e., using the substitution $t=c_1(1)(1-x)$. For $K$ large the function $e^{-Kt}$ goes from the appreciable value $1$ at $t=0$ to the very small value $e^{-\sqrt{K}}$ only a very small distance away at $t=1/\sqrt{K}$. Even more, if one wants the steps in the function values of $e^{-K\cdot n\,dt}$ to be small, one would need $K\,dt$ to be small, $1\gg K\,dt\gg 1/K$, for instance by choosing $dt=K^{-3/2}$.
The mesh refinement algorithm would have to detect this situation. However, without telling the algorithm about this situation (for instance by structuring the initial mesh this way), it would have to create a sizable number, up to $O(K)$, of mesh points with a distance of $K^{-1-\epsilon}$ close to $x=1$ to zero-in to the singularity. For $K=10^5$ this would requite mesh points in the order of some thousands or ten thousands, which will likely exceed any moderate size of Nmax.
It is not appropriate to work in the space $H^1_0(U)$ since nonzero boundary conditions are being considered. You will have to assume some regularity of the boundary of $U$.
One version of Green's theorem (see e.g. the appendix in Evans) is that $$- \int_U (\Delta u) v \, dx = \int_U Du \cdot Dv \, dx -\int_{\partial U} \frac{\partial u}{\partial \nu} v \, dS.$$
A weak solution to the problem at hand can be proposed by setting $-\Delta u = f$ in $U$ and $\dfrac{\partial u}{\partial \nu} = -u$ on $\partial U$ so that $$\int_U fv = \int_U Du \cdot Dv + \int_{\partial U} uv \, dS \quad \forall v \in H^1(U)$$
or a bit more precisely
$$\int_U fv = \int_U Du \cdot Dv + \int_{\partial U} \newcommand{\tr}{\mathrm{Tr}\ \! }( \tr u )(\tr v) \, dS \quad \forall v \in H^1(U)$$
where $\tr : H^1(U) \to L^2(\partial U)$ is the trace operator.
An appropriate bilinear form is thus given by $$B[u,v] = \int_U Du \cdot Dv + \int_{\partial U} \newcommand{\tr}{\mathrm{Tr}\ \! }( \tr u )(\tr v) \, dS, \quad u,v \in H^1(U).$$
$B$ is clearly bounded. As far as coercivity goes, it may be helpful to use the Rellich-Kondrachov theorem. I can follow up with a hint if you like.
It remains to show that there is a constant $\alpha > 0$ with the property that $\|u\|_{H^1}^2 \le \alpha B[u,u]$ for all $u \in H^1(U)$. This can be proven by contradition. Otherwise, for every $n \ge \mathbb N$ there would exist $u_n \in H^1(U)$ with the property that $\|u_n\|^2_{H^1} > n B[u_n,u_n]$. For each $n$ define $v_n = \dfrac{u_n}{\|u_n\|_{H^1}}$. Then $v_n \in H^1(U)$, $\|v_n\|_{H^1} = 1$, and $B[v_n,v_n] < \dfrac 1n$ and all $n$.
Here we can invoke Rellich-Kondrachov. Since the family $\{v_n\}$ is bounded in the $H^1$ norm, there is a subsequence $\{v_{n_k}\}$ that converges to a limit $v \in L^2(U)$. However, since $\|Dv_{n_k}\|_{L^2}^2 < \dfrac{1}{n_k}$ it is also true that $Dv_{n_k} \to 0$ in $L^2$. Thus for any $\phi \in C_0^\infty(U)$ you have $$\int_U v D \phi \, dx = \lim_{k \to \infty} \int_U v_{n_k} D \phi \, dx = - \lim_{k \to \infty} \int_U D v_{n_k} \phi \, dx = 0.$$ This means $v \in H^1(U)$ and $D v = 0$, from which you can conclude $v_{n_k} \to v$ in $H^1(U)$. Since $\|v_{n_k}\|_{H^1} = 1$ for all $k$ it follows that $\|v\|_{H^1} = 1$ as well.
Next, since $\|\tr v_{n_k}\|_{L^2(\partial U)}^2 < \dfrac{1}{n_k}$ and the trace operator is bounded there is a constant $C$ for which
$$ \|\tr v\|_{L^2(\partial \Omega)} \le \|\tr v - \tr v_{n_k}\|_{L^2(\partial \Omega)} + \|\tr v_{n_k}\|_{L^2(\partial \Omega)} < \frac{1}{n_k} + C \|v - v_{n_k}\|_{H^1(U)}.$$
Let $k \to \infty$ to find that $\tr v = 0$ in $L^2(\partial U)$.
Can you prove that if $v \in H^1(U)$, $Dv = 0$, and $\tr v = 0$, then $v = 0$? Once you have established that fact you arrive at a contradiction, since $v$ also satisfies $\|v\|_{H^1} = 1$. It follows that $B$ is in fact coercive.
Best Answer
Your question is equivalent to the question which boundary value problems (ODEs or PDEs) are well posed, in the sense of Hadamard. To this question alone, a major part of PDE research has been devoted for the last century or so, especially in the case of boundary conditions which are considered to be physically meaningful. To quote Wikipedia's lemma on boundary value problems:
A famous example is the question when the Navier-Stokes equations possess a global, smooth solution; it is one of the one million dollar Clay Institute problems. For more thoughts on the issue of uniqueness for this particular problem, see here. For some recent well-posedness results for the Navier-Stokes equations, see here.
Even in the ODE setting, the question of uniqueness related to boundary conditions is a difficult one. Quoting Scholarpedia's lemma on boundary value problems, which focuses exclusively on the ODE context:
I hope this helps!