Newton's method for systems
Many real problems ask you to solve not one equation but several at once: find the values that make a dozen coupled nonlinear relations all vanish together — think of balancing chemical reactions, finding the equilibrium of a circuit, or fitting a physics model. Newton's method extends beautifully to this setting. The single tangent line becomes a tangent PLANE (the best linear approximation of the whole system), and following it to where it hits zero means solving a linear system each step.
Write the system as a vector equation F(x) = 0, where x is a vector of n unknowns and F returns n values. The local linear model is F(x) is approximately F(x_n) + J(x_n)(x - x_n), where J is the JACOBIAN MATRIX — the n-by-n grid of all partial derivatives, J_ij = partial F_i / partial x_j. Setting the model to zero gives the Newton step: solve the linear system J(x_n) * s = -F(x_n) for the step s, then update x_{n+1} = x_n + s. Note we SOLVE the linear system (with an LU factorization) rather than form the inverse J^-1 — forming the inverse is wasteful and less accurate. Each iteration thus costs one Jacobian evaluation plus one linear solve (about (2/3)n^3 flops). Near a simple root (where J is nonsingular) convergence is still quadratic, just as in one dimension.
The same honesty applies as in 1-D, amplified. Newton for systems needs a good starting guess; far from a solution it can diverge or wander, so practical codes add a line search or trust region (damping the step) to make it globally more reliable. The Jacobian can be expensive or unavailable analytically — you may approximate it by finite differences (n extra function evaluations) or by automatic differentiation, and quasi-Newton methods like Broyden's avoid recomputing it every step. If J is singular or ill-conditioned at the root, convergence degrades exactly as a near-zero derivative hurts the 1-D method.
Solve the two equations x^2 + y^2 = 1 and y = x^2 (a circle meeting a parabola). Write F = [x^2 + y^2 - 1, y - x^2]; the Jacobian is J with rows (2x, 2y) and (-2x, 1). From a guess like (0.8, 0.6), each step solves the 2-by-2 system J s = -F and updates, converging quadratically to the intersection (0.786, 0.618).
Solve J s = -F each step; the Jacobian is the tangent plane of the whole system.
Each step SOLVES J s = -F (an LU factorization), never forms J^-1 — inverting is wasteful and less accurate. And quadratic convergence still demands a good initial guess and a nonsingular Jacobian; far away you need a line search or trust region.