I'm looking for the vector valued solution $y(t)=(y_1(t),\ldots,y_n(t))$ of a linear system of ODEs in the form of
$$y'(t)=f(t,y(t)),\ \ \ \mathrm{where}\ f(t,y)=M\cdot y+g(t)$$
The matrix $M \in \mathbb R^{n \times n}$ and the continuous, vector valued function $g$ are given.
Furthermore we have an initial value, given as a vector:
$$y(0)=y_0, \ \ \ y_0 \in \mathbb R^n$$
Method:
One way to approximate the solution at $t=h$ numerically, is to choose an (appropriate) matrix $A=(a_{ij}) \in \mathbb R^{s \times s}$, a parameter vector $c=(c_i) \in \mathbb R^s$ $(0 \leq c_i \leq 1)$ and a vector $b=(b_i) \in \mathbb R^s$. And write down the following system of $s \cdot n$ equations and $s \cdot n$ unknowns:
\begin{align}\tag{1} &y_1=y_0+h(a_{11}f(c_1h,y_1)+\ldots +a_{1s}f(c_sh,y_s))\\&\vdots\\&y_s=y_0+h(a_{s1}f(c_1h,y_1)+\ldots +a_{ss}f(c_sh,y_s))\end{align}
where $y_i=(y_{i,1},\ldots,y_{i,n}) \in \mathbb R^n, i=1,\ldots,s$ have to be determined. After this system is solved, an approximation of the solution $y(t)$ at $t=h$ is given by
$$\tag{2}y_0+h(b_1f(c_1h,y_1)+\ldots+b_sf(c_sh,y_s)$$
How can I show that the system of equations (1) has an unique solution for sufficiently small stepsize $h$?
Furthermore, does this method have a specific name? I'm also looking for a matlab code, that solves (1) and evaluates (2)
I don't understand your definition of $F$. Don't we have $Y=Y_0+hAF(T,Y)$ with $F(T,Y)$ defined as
$F(T,Y):=\begin{pmatrix}f(c_1h,y_1) & \dots & f(c_1h,y_1)\\ \vdots & \ddots &\vdots \\ f(c_sh,y_s) & \dots & f(c_sh,y_s) \end{pmatrix}$.
Furthermore, why is $Y=Y_0+hAF(T,Y)$ a fixed-point equation? What exactly does fixed point equation mean? In Banachs theorem we have $x_{n+1}=\phi(x_n)$.
Do we need $h$ small enough such that $hL\max_{i}\sum_j|a_{ij}|<1$? (So that we have a contracting map)
Best Answer
General argument based on the implicit function theorem
The implicit system that you have to solve is $$ y_i=y_0+h\bigl(a_{i1}f(t_1,y-1)+...+a_{is}f(t_s,y_s)\bigr).\tag{1a} $$ This can be written as one big system $$Y=Y_0+h(A\otimes I)F(T,Y),\tag{1b}$$ where $$ Y=\pmatrix{y_1\\\vdots\\y_s}, ~~Y_0=\pmatrix{y_0\\\vdots\\y_0}, ~~ F(T,Y)=\pmatrix{f(t_1,y_1)\\\vdots\\f(t_s,y_s)} $$ and the tensor product is the Kronecker product. For $h=0$ this has the obvious and unique solution $Y=Y_0$. The composite map $F$ is as smooth as the given function $f$, and $Y\mapsto Y-h(A\otimes I)F(T,Y)$ has an invertible Jacobian at $Y=Y_0$ for sufficiently small $h$. The existence and uniqueness of the solution of equation (1) in dependence of $h$ in some small interval around $0$ can now be asserted via the implicit function theorem.
Direct solution in the linear case
If $f(t,y)=My+g(t)$ as per question, then the implicit system is in fact linear and can be solved as an $ns\times ns$ system, $$ Y=Y_0+h(A\otimes I)[(I\otimes M)Y+G(T)] \\\implies\\ Y=[I-h(A\otimes M)]^{-1}[Y_0+h(A\otimes I)G(T)] $$ where $G(T)$ is again the vertically stacked collection of vectors $(g(t_1),...,g(t_s))$.
As no further numerical method is involved, the only restriction on $h$ is that $h^{-1}$ should not be an eigenvalue of $(A\otimes M)$. This is satisfied if $h\|A\|\,\|M\|<1$.
Using the Banach fixed-point theorem directly for the non-linear case
One can make this a bit more quantitative using the Banach fixed-point theorem for the fixed-point map $$P(Y)=Y_0+h(A\otimes I)F(T,Y).$$ Using a suitable composite norm on $Y$, for instance $\|Y\|=\max_{k=1,..,s}\|y_k\|$, one finds that $P$ has a Lipschitz constant $$L(P)=h\|A\|\,L(f)$$ on some ball $\|Y-Y_0\|\le R$ where $R$ is such that $f$ has Lipschitz constant $L(f)$ on $\|y-y_0\|\le R$. Choosing $h$ small enough, one can achieve $L(P)\le\frac12$.
Then with $$\|P(Y)-Y_0\|\le \|P(Y_0)-Y_0\|+\|P(Y)-P(Y_0)\|\le h\|A\|\|F(T,Y_0)\|+L(P)\|Y-Y_0\|$$ one needs to ensure that also the first term is smaller than $R/2$ to have $\|Y-Y_0\|\le R\implies \|P(Y)-Y_0\|\le R$ so that all conditions of the fixed-point theorem are satisfied.
Some classes of implicit methods
Examples that you might want to explore are the Gauß methods and the SDIRK methods. The Gauß methods, based on Gauß integration, have the highest order for the number of stages. The SDIRK methods (Singly Diagonally Implicit Runge-Kutta) have a lower triangular matrix $(a_{ij})$ with the same entry on the diagonal. This has the advantage that the implicit system can be solved stage-by-stage, and in each stage a simplified Newton method can be employed that has the same matrix.
Butcher tableaus for some lower order examples can be found in these slides of Butcher's.