# Solving system of nonlinear vector functions

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.

• For the Newton method you'll need write your equations as $\vec{F}(\vec{x})=0$ where $\vec{x}$ is the vector of unknowns. Then this will be the familiar situation of a single vector function. – Maxim Umansky Sep 5 '20 at 16:35
• @MaximUmansky I'm aware of that, but I am still confused on how to properly set up a system in $\mathbb{R}^{2m}$ – vlovero Sep 5 '20 at 16:43

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.