In order to solve a boundary problem for a conducting hemisphere, the following matrix equation arises (derived from the boundary condition on the curved part of the surface) where we must solve for the quantities $b_n$, $n\in\Bbb N$: \begin{align} \sum_{n=0}^\infty U_{nk} b_n = a_k \label{1}\tag{1} \end{align} for $k=0,1,2,...\infty$, where $a_k$ are known and the matrix elements $U_{nk}$ are $$\begin{cases} U_{n,n}&=\frac{1}{2n+1} \\ U_{2n,2k}&=0 \qquad (n\neq k)\\ U_{2n+1,2k}&=-\frac{(2n+1)P_{2n}(0)P_{2k}(0)}{2(2k-2n-1)(n+k+1)} \\ U_{2n,2k+1}&=-\frac{(2k+1)P_{2n}(0)P_{2k}(0)}{2(2n-2k-1)(n+k+1)} \end{cases} $$ and $P_n(0)=(-1)^{n/2}(n-1)!!/n!!$ are the Legendre polynomials evaluated at 0. Note the matrix for $U_{nk}$ is symmetric, and checkered - all entries of the same parity are zero except the diagonal. (I don't know how relevant the exact values of the matrix elements are to my question of solveability, but there they are anyway).
I know that Eq. \eqref{1} is ill-posed/underdetermined because it only uses one of the boundary conditions; the true solution for $b_n$ can be found by combining \eqref{1} with other equations to obtain a final matrix equation, which I will name (2). What confuses me is that we can still solve Eq. \eqref{1} by itself, and apparently get a unique answer. For example if $a_k=U_{0k}$, then simply $b_n=\delta_{n0}$ solves \eqref{1}. But the solution to (2) also solves \eqref{1}, and I imagine that the other boundary conditions could be many different things, all resulting in different equations for (2), and hence different solutions for \eqref{1}.
My question is: how can we tell that Eq. \eqref{1} is underdetermined?
Also why might Matlab simply give one solution for \eqref{1} without a warning of non-uniqueness?