The question
Consider a linear system of size $10 \times 10$ in the form $\mathbf{D}(\alpha) \underline{R}(\alpha) = \underline{L}(\alpha)$, with $\alpha \in [0,A]$ and $A$ being a large real number. The entries of $D_{ij}(\alpha) = D(\alpha,m_j(\alpha))$ are complicated complex-valued functions of the complex polynomial roots $m_j(\alpha)$, which I cannot derive in closed form but I have computed numerically. From inspection and from numerical evaluation, I know that:
$$ \Re(D_{i,j}) = - \Re(D_{i,j+4}); \;\; \Im(D_{i,j}) = \Im(D_{i,j+4})$$ for $i = 1,2,4,5,8; j \leq 4$, and $$ \Re(D_{i,j}) = \Re(D_{i,j+4}); \;\; \Im(D_{i,j}) = - \Im(D_{i,j+4})$$
for $i = 3,6,9; j \leq 4$. For $I=7,10$ the coefficients are real constants. Does this structure, i.e. reflection around the imaginary/real axes, somehow entail ill-conditioning? My guess is that this property could be seen as a linear dependence of the columns, and therefore it could be the source of the matrix condition number having $\min{\kappa(\alpha)} \approx$ 1e+4 and, more worryingly, $\max{\kappa(\alpha)} \approx$ 1e+12. I performed an asymptotic analysis for $\alpha \to 0$ and $\mathbf{D}(\alpha)$ becomes singular, while for large $A$(i.e., $\alpha \to \infty$) I know the matrix becomes ill conditioned from numerical calculations.
If this is the case, are there known approaches to improve the conditioning and obtain more precise solutions? I have resorted to using quadruple precision in Matlab, but I would prefer a better analytical justification.
Details for the interested reader
The problem arises from a complicated boundary value in the context of fracture mechanics of thin shells. I am solving a set of coupled PDEs with mixed boundary conditions, to which I apply Fourier transforms to obtain a set of ODEs. Assuming the transform to have the form of a superposition of exponentials $f(x,\alpha) = R_j(\alpha) e^{m_j x}$, the problem can be reduced to the solution of an eigenvalue problem. Given the expected singularity in the result, I apply a variable transformation to the problem unknowns $R_j = \sum_k B_{ik} q_k$, where $B_{ik}$ is the result of $\mathbf{D}(\alpha) (B_{ik} q_k) = \underline{L}(\alpha,q_k)$.
This transformation requires the mapping of the boundary conditions to the Fourier space, which consists in the solution of the linear system I mentioned. I asked a related question in the computational stack exchange here.