I have an ODE such as the one below,
$$y"=a^{2}y$$
The analytical solution is clearly
$$y=C_{1}e^{ax} + C_{2}e^{-ax}$$
To apply boundary conditions $y(x=x_{1})=A$ and $y(x=x_{2})=B$, I recast into a matrix:
$\left( \begin{array}{ccc}e^{ax_{1}} & e^{-ax_{1}} \\ e^{ax_{2}} & e^{-ax_{2}} \end{array} \right)$$\left( \begin{array}{ccc}C_{1}\\ C_{2}\end{array} \right)=\left( \begin{array}{ccc}A\\ B\end{array} \right)$
My issue is that when $x_{1}$ and $x_{2}$ are far apart from one another (say $a=1,x_{1}=-10^{4},x_{2}=10^{4}$ or when $a$ is large (say $a=10^{4},x_{1}=-1,x_{2}=1$), the matrix inversion cannot be solved within machine precision. The issue is the exponential, which varies so widely. I was wondering what steps I could take to find the constants to satisfy my ODE solution.
Perhaps it might help to write the solution as
$$ y(x) = \frac{\sinh(a (x-x_2))}{\sinh(a (x_1 - x_2))} A - \frac{\sinh(a (x - x_1))}{\sinh(a (x_1 - x_2))} B$$
There will still be problems with overflow/underflow in numerical computation, though. Of course, when $t$ is large and positive, $\sinh(t) \sim e^t/2$, while when $t$ is large and negative, $\sinh(t) \sim -e^{-t}/2$. These can be used to get useful approximations in various regions.