Solving system of nonlinear vector functions

Question

I am trying to figure out how to implement a solver for a system of nonlinear equations of the form

\begin{align*} u_1 &= y_n + h\left(a_{1,1}f(t_n + c_1 h, u_1) + a_{1,2}f(t_n + c_2 h, u_2)\right) \\ u_2 &= y_n + h\left(a_{2,1}f(t_n + c_1 h, u_1) + a_{2,2}f(t_n + c_2 h, u_2)\right) \end{align*} where $f: \mathbb{R}^m\to\mathbb{R}^m$ can be any nonlinear function, and $u_1, u_2 \in \mathbb{R}^m$ are the only unknowns.

I know how to use Newton's method $\vec{x}_{k+1} = \vec{x}_{k} - J^{-1}(F)F(\vec{x}_{k})$ for a single vector function, but I am confused on how to adapt this for multiple. From the papers I have been reading, authors reference using a modified Newton method for block matrices created with the kronecker product, but when I do that it leaves me with a matrix in $\mathbb{R}^{2m}$ that I don't know what to do with. I have also seen authors define a Jacobian matrix that contains other Jacobians, but again, I don't know how to handle that on a computer. How should I go about creating an iterative method for a system like this? I am trying to implement this to use it with the Radau IIA methods.

Answer 1

If you have a single vector equation $\vec{F}(\vec{x})=0$ then you solve it by representing that state vector $\vec{x}$ as a set of amplitudes $[x_0,x_1,...,x_{n-1}]$ after discretization by your favorite method (FD, FV, FEM, spectral); and we know how to solve it. If you also have a second equation $\vec{G}(\vec{y})=0$ then the full state vector is $[x_0,x_1,...,x_{n-1}, y_0,y_1, ..., y_{n-1}]$, otherwise the same thing.