# Efficiently solving a large, sparse linear system $M(s)ab(s)=c(s)$ (determined by smooth functions) over some range of $s$

I'm looking at a differential equation on the edges of a graph (the application is neuroscience), and the Laplace transform of the solution on most of the edges has a general solution more-or-less of the form $(\mathcal{L}v)(x,s)= a(s)e^{-x/\sqrt{s}} + b(s)e^{x/\sqrt{s}}$, where $x$ is in the interval assigned to the edge and $s$ is the Laplace variable. One of the edges has another term in its solution which has no unknowns and can be evaluated easily.

The functions $a(s)$ and $b(s)$ are unknown, but can be found by putting together the boundary conditions for each vertex of my graph. This means substituting $x=0$ or $x=L$ (depending on the direction of the edge) for each edge incident on each vertex, and substituting a value of $s$ ($L$ is the length of an edge). This results in a system of $2m$ independent linear equations in $2m$ variables (where $m$ is the number of edges). We can write this system as $M(s)ab(s)=c(s)$, where $M(s)$ is a sparse matrix whose nonzero entries are smooth functions of $s$, $c(s)$ is a vector mostly consisting of zeros whose entries are smooth functions of $s$, and the unknown $ab(s)$ contains the values of $a(s)$ and $b(s)$ for all edges. The zero entries are zero for every value of $s$.

I plan to get the solution $v(x,t)$ for a single value of $x$ on an edge (in the time domain) by finding $a(s)$ and $b(s)$ for that edge, for some range of $s$ (probably a large finite number of discrete values), substituting $a(s)$ and $b(s)$ into the formula for $(\mathcal{L}v)(x,s)$, and feeding the resulting Laplace-domain solution into a numerical inverse Laplace transform (using a fast Fourier transform) to get a time series of $v(x,t)$ (for a single value of $x$ and a range of $t$-values).

An obvious way to do this is to apply a sparse linear system solver to the system $S$ times, where $S$ is the number of values of $s$ required for the NILT. This could be very time-consuming, and I would like to be able to do this for values of $m$ as large as possible. My question(s) are:

Is there a name for this kind of problem? I would be surprised if no-one has encountered anything similar before, but I've no idea what it would be called. EDIT: These are called parametrised matrix equations. Most of the material about them that I can find is by Paul Constantine.

Can the smoothness of the functions defining $M(s)$ and $c(s)$ be exploited to find the array of values of $a(s)$ and $b(s)$ without solving the whole system $S$ times? For instance, is it possible to specify any the circumstances under which interpolation could be used to get intermediate values of $a(s)$ and $b(s)$ without ruining the accuracy of the NILT at the end?

Would it be practical to solve this symbolically using a computer algebra system? All of the nonzero matrix entries have a form something like $Ke^{L/\sqrt{s}}$, where $K$ and $L$ are known numbers (and $L$ may be zero).