View source for Multidimensional Newton-Raphson method

{{Draft task|Arithmetic operations}}

;Task:
Create a program that finds and outputs the root of a system of nonlinear equations
using Newton-Raphson method.
<br><br>

=={{header|C sharp|C#}}==
For matrix inversion and matrix and vector definitions - see C# source from [[Gaussian elimination]]
<syntaxhighlight lang="csharp">
using System;

namespace Rosetta
{
    internal interface IFun
    {
        double F(int index, Vector x);
        double df(int index, int derivative, Vector x);
        double[] weights();
    }

    class Newton
    {                
        internal Vector Do(int size, IFun fun, Vector start)
        {
            Vector X = start.Clone();
            Vector F = new Vector(size);
            Matrix J = new Matrix(size, size);
            Vector D;
            do
            {
                for (int i = 0; i < size; i++)
                    F[i] = fun.F(i, X);
                for (int i = 0; i < size; i++)
                    for (int j = 0; j < size; j++)
                        J[i, j] = fun.df(i, j, X);
                J.ElimPartial(F);
                X -= F;
                //need weight vector because different coordinates can diffs by order of magnitudes
            } while (F.norm(fun.weights()) > 1e-12);
            return X;
        }       
    }
}
</syntaxhighlight>
<syntaxhighlight lang="csharp">
using System;

//example from https://eti.pg.edu.pl/documents/176593/26763380/Wykl_AlgorOblicz_7.pdf
namespace Rosetta
{
    class Program
    {
        class Fun: IFun
        {
            private double[] w = new double[] { 1,1 };

            public double F(int index, Vector x)
            {
                switch (index)
                {
                    case 0: return Math.Atan(x[0]) - x[1] * x[1] * x[1];
                    case 1: return 4 * x[0] * x[0] + 9 * x[1] * x[1] - 36;
                }
                throw new Exception("bad index");
            }

            public double df(int index, int derivative, Vector x)
            {
                switch (index)
                {
                    case 0:
                        switch (derivative)
                        {
                            case 0: return 1 / (1 + x[0] * x[0]);
                            case 1: return -3*x[1]*x[1];
                        }
                        break;
                    case 1:
                        switch (derivative)
                        {
                            case 0: return 8 * x[0];
                            case 1: return 18 * x[1];
                        }
                        break;
                }
                throw new Exception("bad index");
            }
            public double[] weights() { return w; }
        }

        static void Main(string[] args)
        {
            Fun fun = new Fun();
            Newton newton = new Newton();
            Vector start = new Vector(new double[] { 2.75, 1.25 });
            Vector X = newton.Do(2, fun, start);
            X.print();
        }
    }
}
</syntaxhighlight>
{{out}}<pre>
2.54258545959024
1.06149981539336
</pre>

=={{header|Go}}==
{{trans|Kotlin}}
<br>
We follow the Kotlin example of coding our own matrix methods rather than using a third party library.
<syntaxhighlight lang="go">package main

import (
    "fmt"
    "math"
)

type vector = []float64
type matrix []vector
type fun = func(vector) float64
type funs = []fun
type jacobian = []funs

func (m1 matrix) mul(m2 matrix) matrix {
    rows1, cols1 := len(m1), len(m1[0])
    rows2, cols2 := len(m2), len(m2[0])
    if cols1 != rows2 {
        panic("Matrices cannot be multiplied.")
    }
    result := make(matrix, rows1)
    for i := 0; i < rows1; i++ {
        result[i] = make(vector, cols2)
        for j := 0; j < cols2; j++ {
            for k := 0; k < rows2; k++ {
                result[i][j] += m1[i][k] * m2[k][j]
            }
        }
    }
    return result
}

func (m1 matrix) sub(m2 matrix) matrix {
    rows, cols := len(m1), len(m1[0])
    if rows != len(m2) || cols != len(m2[0]) {
        panic("Matrices cannot be subtracted.")
    }
    result := make(matrix, rows)
    for i := 0; i < rows; i++ {
        result[i] = make(vector, cols)
        for j := 0; j < cols; j++ {
            result[i][j] = m1[i][j] - m2[i][j]
        }
    }
    return result
}

func (m matrix) transpose() matrix {
    rows, cols := len(m), len(m[0])
    trans := make(matrix, cols)
    for i := 0; i < cols; i++ {
        trans[i] = make(vector, rows)
        for j := 0; j < rows; j++ {
            trans[i][j] = m[j][i]
        }
    }
    return trans
}

func (m matrix) inverse() matrix {
    le := len(m)
    for _, v := range m {
        if len(v) != le {
            panic("Not a square matrix")
        }
    }
    aug := make(matrix, le)
    for i := 0; i < le; i++ {
        aug[i] = make(vector, 2*le)
        copy(aug[i], m[i])
        // augment by identity matrix to right
        aug[i][i+le] = 1
    }
    aug.toReducedRowEchelonForm()
    inv := make(matrix, le)
    // remove identity matrix to left
    for i := 0; i < le; i++ {
        inv[i] = make(vector, le)
        copy(inv[i], aug[i][le:])
    }
    return inv
}

// note: this mutates the matrix in place
func (m matrix) toReducedRowEchelonForm() {
    lead := 0
    rowCount, colCount := len(m), len(m[0])
    for r := 0; r < rowCount; r++ {
        if colCount <= lead {
            return
        }
        i := r

        for m[i][lead] == 0 {
            i++
            if rowCount == i {
                i = r
                lead++
                if colCount == lead {
                    return
                }
            }
        }

        m[i], m[r] = m[r], m[i]
        if div := m[r][lead]; div != 0 {
            for j := 0; j < colCount; j++ {
                m[r][j] /= div
            }
        }

        for k := 0; k < rowCount; k++ {
            if k != r {
                mult := m[k][lead]
                for j := 0; j < colCount; j++ {
                    m[k][j] -= m[r][j] * mult
                }
            }
        }
        lead++
    }
}

func solve(fs funs, jacob jacobian, guesses vector) vector {
    size := len(fs)
    var gu1 vector
    gu2 := make(vector, len(guesses))
    copy(gu2, guesses)
    jac := make(matrix, size)
    for i := 0; i < size; i++ {
        jac[i] = make(vector, size)
    }
    tol := 1e-8
    maxIter := 12
    iter := 0
    for {
        gu1 = gu2
        g := matrix{gu1}.transpose()
        t := make(vector, size)
        for i := 0; i < size; i++ {
            t[i] = fs[i](gu1)
        }
        f := matrix{t}.transpose()
        for i := 0; i < size; i++ {
            for j := 0; j < size; j++ {
                jac[i][j] = jacob[i][j](gu1)
            }
        }
        g1 := g.sub(jac.inverse().mul(f))
        gu2 = make(vector, size)
        for i := 0; i < size; i++ {
            gu2[i] = g1[i][0]
        }
        iter++
        any := false
        for i, v := range gu2 {
            if math.Abs(v)-gu1[i] > tol {
                any = true
                break
            }
        }
        if !any || iter >= maxIter {
            break
        }
    }
    return gu2
}

func main() {
    /*
       solve the two non-linear equations:
       y = -x^2 + x + 0.5
       y + 5xy = x^2
       given initial guesses of x = y = 1.2

       Example taken from:
       http://www.fixoncloud.com/Home/LoginValidate/OneProblemComplete_Detailed.php?problemid=286

       Expected results: x = 1.23332, y = 0.2122
    */
    f1 := func(x vector) float64 { return -x[0]*x[0] + x[0] + 0.5 - x[1] }
    f2 := func(x vector) float64 { return x[1] + 5*x[0]*x[1] - x[0]*x[0] }
    fs := funs{f1, f2}
    jacob := jacobian{
        funs{
            func(x vector) float64 { return -2*x[0] + 1 },
            func(x vector) float64 { return -1 },
        },
        funs{
            func(x vector) float64 { return 5*x[1] - 2*x[0] },
            func(x vector) float64 { return 1 + 5*x[0] },
        },
    }
    guesses := vector{1.2, 1.2}
    sol := solve(fs, jacob, guesses)
    fmt.Printf("Approximate solutions are x = %.7f,  y = %.7f\n", sol[0], sol[1])

    /*
       solve the three non-linear equations:
       9x^2 + 36y^2 + 4z^2 - 36 = 0
       x^2 - 2y^2 - 20z = 0
       x^2 - y^2 + z^2 = 0
       given initial guesses of x = y = 1.0 and z = 0.0

       Example taken from:
       http://mathfaculty.fullerton.edu/mathews/n2003/FixPointNewtonMod.html (exercise 5)

       Expected results: x = 0.893628, y = 0.894527, z = -0.0400893
    */

    fmt.Println()
    f3 := func(x vector) float64 { return 9*x[0]*x[0] + 36*x[1]*x[1] + 4*x[2]*x[2] - 36 }
    f4 := func(x vector) float64 { return x[0]*x[0] - 2*x[1]*x[1] - 20*x[2] }
    f5 := func(x vector) float64 { return x[0]*x[0] - x[1]*x[1] + x[2]*x[2] }
    fs = funs{f3, f4, f5}
    jacob = jacobian{
        funs{
            func(x vector) float64 { return 18 * x[0] },
            func(x vector) float64 { return 72 * x[1] },
            func(x vector) float64 { return 8 * x[2] },
        },
        funs{
            func(x vector) float64 { return 2 * x[0] },
            func(x vector) float64 { return -4 * x[1] },
            func(x vector) float64 { return -20 },
        },
        funs{
            func(x vector) float64 { return 2 * x[0] },
            func(x vector) float64 { return -2 * x[1] },
            func(x vector) float64 { return 2 * x[2] },
        },
    }
    guesses = vector{1, 1, 0}
    sol = solve(fs, jacob, guesses)
    fmt.Printf("Approximate solutions are x = %.7f,  y = %.7f,  z = %.7f\n", sol[0], sol[1], sol[2])
}</syntaxhighlight>

{{out}}
<pre>
Approximate solutions are x = 1.2333178,  y = 0.2122450

Approximate solutions are x = 0.8936282,  y = 0.8945270,  z = -0.0400893
</pre>

=={{header|Julia}}==
NLsolve is a Julia package for nonlinear systems of equations, with the Newton-Raphson method one of the choices for solvers.
<syntaxhighlight lang="julia"># from the NLSolve documentation: to solve
#     (x, y) -> ((x+3)*(y^3-7)+18, sin(y*exp(x)-1))
using NLsolve

function f!(F, x)
    F[1] = (x[1]+3)*(x[2]^3-7)+18
    F[2] = sin(x[2]*exp(x[1])-1)
end

function j!(J, x)
    J[1, 1] = x[2]^3-7
    J[1, 2] = 3*x[2]^2*(x[1]+3)
    u = exp(x[1])*cos(x[2]*exp(x[1])-1)
    J[2, 1] = x[2]*u
    J[2, 2] = u
end

println(nlsolve(f!, j!, [ 0.1; 1.2], method = :newton))
</syntaxhighlight>{{out}}
<pre>
Results of Nonlinear Solver Algorithm
 * Algorithm: Newton with line-search
 * Starting Point: [0.1, 1.2]
 * Zero: [-3.7818e-16, 1.0]
 * Inf-norm of residuals: 0.000000
 * Iterations: 4
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 5
 * Jacobian Calls (df/dx): 4
</pre>

=={{header|Kotlin}}==
A straightforward approach multiplying by the inverse of the Jacobian, rather than dividing by f'(x) as one would do in the single dimensional case, which is quick enough here.

As neither the JDK nor the Kotlin Standard Library have matrix functions built in, most of the functions used have been taken from other tasks. 
<syntaxhighlight lang="scala">// Version 1.2.31

import kotlin.math.abs

typealias Vector = DoubleArray
typealias Matrix = Array<Vector>
typealias Func = (Vector) -> Double
typealias Funcs = List<Func>
typealias Jacobian = List<Funcs>

operator fun Matrix.times(other: Matrix): Matrix {
    val rows1 = this.size
    val cols1 = this[0].size
    val rows2 = other.size
    val cols2 = other[0].size
    require(cols1 == rows2)
    val result = Matrix(rows1) { Vector(cols2) }
    for (i in 0 until rows1) {
        for (j in 0 until cols2) {
            for (k in 0 until rows2) {
                result[i][j] += this[i][k] * other[k][j]
            }
        }
    }
    return result
}

operator fun Matrix.minus(other: Matrix): Matrix {
    val rows = this.size
    val cols = this[0].size
    require(rows == other.size && cols == other[0].size)
    val result = Matrix(rows) { Vector(cols) }
    for (i in 0 until rows) {
        for (j in 0 until cols) {
            result[i][j] = this[i][j] - other[i][j]
        }
    }
    return result
} 

fun Matrix.transpose(): Matrix {
    val rows = this.size
    val cols = this[0].size
    val trans = Matrix(cols) { Vector(rows) }
    for (i in 0 until cols) {
        for (j in 0 until rows) trans[i][j] = this[j][i]
    }
    return trans
}

fun Matrix.inverse(): Matrix {
    val len = this.size
    require(this.all { it.size == len }) { "Not a square matrix" }
    val aug = Array(len) { DoubleArray(2 * len) }
    for (i in 0 until len) {
        for (j in 0 until len) aug[i][j] = this[i][j]
        // augment by identity matrix to right
        aug[i][i + len] = 1.0
    }
    aug.toReducedRowEchelonForm()
    val inv = Array(len) { DoubleArray(len) }
    // remove identity matrix to left
    for (i in 0 until len) {
        for (j in len until 2 * len) inv[i][j - len] = aug[i][j]
    }
    return inv
}

fun Matrix.toReducedRowEchelonForm() {
    var lead = 0
    val rowCount = this.size
    val colCount = this[0].size
    for (r in 0 until rowCount) {
        if (colCount <= lead) return
        var i = r

        while (this[i][lead] == 0.0) {
            i++
            if (rowCount == i) {
                i = r
                lead++
                if (colCount == lead) return
            }
        }

        val temp = this[i]
        this[i] = this[r]
        this[r] = temp

        if (this[r][lead] != 0.0) {
           val div = this[r][lead]
           for (j in 0 until colCount) this[r][j] /= div
        }

        for (k in 0 until rowCount) {
            if (k != r) {
                val mult = this[k][lead]
                for (j in 0 until colCount) this[k][j] -= this[r][j] * mult
            }
        }

        lead++
    }
}

fun solve(funcs: Funcs, jacobian: Jacobian, guesses: Vector): Vector {
    val size = funcs.size
    var gu1: Vector
    var gu2 = guesses.copyOf()
    val jac = Matrix(size) { Vector(size) }
    val tol = 1.0e-8
    val maxIter = 12
    var iter = 0
    do {
        gu1 = gu2
        val g = arrayOf(gu1).transpose()
        val f = arrayOf(Vector(size) { funcs[it](gu1) }).transpose()
        for (i in 0 until size) {
            for (j in 0 until size) {
                jac[i][j] = jacobian[i][j](gu1)
            }
        }
        val g1 = g - jac.inverse() * f
        gu2 = Vector(size) { g1[it][0] }
        iter++
    }
    while (gu2.withIndex().any { iv -> abs(iv.value - gu1[iv.index]) > tol } && iter < maxIter)
    return gu2
}

fun main(args: Array<String>) {
    /* solve the two non-linear equations:
       y = -x^2 + x + 0.5
       y + 5xy = x^2
       given initial guesses of x = y = 1.2

       Example taken from:
       http://www.fixoncloud.com/Home/LoginValidate/OneProblemComplete_Detailed.php?problemid=286

       Expected results: x = 1.23332, y = 0.2122
    */

    val f1: Func = { x -> -x[0] * x[0] + x[0] + 0.5 - x[1] }
    val f2: Func = { x -> x[1] + 5 * x[0] * x[1] - x[0] * x[0] }
    val funcs = listOf(f1, f2)
    val jacobian = listOf(
        listOf<Func>({ x -> - 2.0 * x[0] + 1.0 }, { _ -> -1.0 }),
        listOf<Func>({ x -> 5.0 * x[1] - 2.0 * x[0] }, { x -> 1.0 + 5.0 * x[0] })
    )
    val guesses = doubleArrayOf(1.2, 1.2)
    val (xx, yy) = solve(funcs, jacobian, guesses)
    System.out.printf("Approximate solutions are x = %.7f,  y = %.7f\n", xx, yy)

    /* solve the three non-linear equations:
       9x^2 + 36y^2 + 4z^2 - 36 = 0
       x^2 - 2y^2 - 20z = 0
       x^2 - y^2 + z^2 = 0
       given initial guesses of x = y = 1.0 and z = 0.0

       Example taken from:
       http://mathfaculty.fullerton.edu/mathews/n2003/FixPointNewtonMod.html (exercise 5)

       Expected results: x = 0.893628, y = 0.894527, z = -0.0400893
    */

    println()
    val f3: Func = { x -> 9.0 * x[0] * x[0] + 36.0 * x[1] * x[1] + 4.0 * x[2] * x[2] - 36.0 }
    val f4: Func = { x -> x[0] * x[0] - 2.0 * x[1] * x[1] - 20.0 * x[2] }
    val f5: Func = { x -> x[0] * x[0] - x[1] * x[1] + x[2] * x[2] }
    val funcs2 = listOf(f3, f4, f5)
    val jacobian2 = listOf(
        listOf<Func>({ x -> 18.0 * x[0] }, { x -> 72.0 * x[1] }, { x -> 8.0 * x[2] }),
        listOf<Func>({ x -> 2.0 * x[0] }, { x -> -4.0 * x[1] }, { _ -> -20.0 }),
        listOf<Func>({ x -> 2.0 * x[0] }, { x -> -2.0 * x[1] }, { x -> 2.0 * x[2] })
    )
    val guesses2 = doubleArrayOf(1.0, 1.0, 0.0)
    val (xx2, yy2, zz2) = solve(funcs2, jacobian2, guesses2)
    System.out.printf("Approximate solutions are x = %.7f,  y = %.7f,  z = %.7f\n", xx2, yy2, zz2)
}</syntaxhighlight>

{{out}}
<pre>
Approximate solutions are x = 1.2333178,  y = 0.2122450

Approximate solutions are x = 0.8936282,  y = 0.8945270,  z = -0.0400893
</pre>

=={{header|Nim}}==
{{trans|Kotlin}}
<syntaxhighlight lang="nim">import sequtils, strformat, sugar

type
  Vector = seq[float]
  Matrix = seq[Vector]
  Func = (Vector) -> float
  Funcs = seq[Func]
  Jacobian = seq[Funcs]


func `*`(m1, m2: Matrix): Matrix =
  let
    rows1 = m1.len
    cols1 = m1[0].len
    rows2 = m2.len
    cols2 = m2[0].len
  doAssert cols1 == rows2
  result = newSeqWith(rows1, newSeq[float](cols2))
  for i in 0..<rows1:
    for j in 0..<cols2:
      for k in 0..<rows2:
        result[i][j] += m1[i][k] * m2[k][j]


func `-`(m1, m2: Matrix): Matrix =
  let
    rows = m1.len
    cols = m1[0].len
  doAssert m2.len == rows and m2[0].len == cols
  result = newSeqWith(rows, newSeq[float](cols))
  for i in 0..<rows:
    for j in 0..<cols:
      result[i][j] = m1[i][j] - m2[i][j]


func transposed(m: Matrix): Matrix =
  let
    rows = m.len
    cols = m[0].len
  result = newSeqWith(cols, newSeq[float](rows))
  for i in 0..<cols:
    for j in 0..<rows:
      result[i][j] = m[j][i]


func toReducedRowEchelonForm(m: var Matrix) =
  var lead = 0
  let rowCount = m.len
  let colCount = m[0].len
  for r in 0..<rowCount:
    if colCount <= lead: return
    var i = r

    while m[i][lead] == 0:
      inc i
      if rowCount == i:
        i = r
        inc lead
        if colCount == lead: return

    swap m[i], m[r]

    if m[r][lead] != 0:
      let divisor = m[r][lead]
      for j in 0..<colCount: m[r][j] /= divisor

    for k in 0..<rowCount:
      if k != r:
        let mult = m[k][lead]
        for j in 0..<colCount: m[k][j] -= m[r][j] * mult

    inc lead


func inverse(m: Matrix): Matrix =
  let size = m.len
  doAssert m.allIt(it.len == size), "not a square matrix."
  var aug = newSeqWith(size, newSeq[float](2 * size))
  for i in 0..<size:
    for j in 0..<size: aug[i][j] = m[i][j]
    # Augment by identity matrix to right.
    aug[i][i + size] = 1
  aug.toReducedRowEchelonForm()
  result = newSeqWith(size, newSeq[float](size))
  # Remove identity matrix to left.
  for i in 0..<size:
    for j in 0..<size: result[i][j] = aug[i][j + size]


proc solve(funcs: Funcs; jacobian: Jacobian; guesses: Vector): Vector =
  let size = funcs.len
  result = guesses
  var jac = newSeqWith(size, newSeq[float](size))
  const Tol = 1e-8
  let MaxIter = 12
  var iter = 1
  while true:
    let gu = move(result)
    let g = transposed(@[gu])
    let f = transposed(@[funcs.mapIt(it(gu))])
    for i in 0..<size:
      for j in 0..<size:
        jac[i][j] = jacobian[i][j](gu)
    let g1 = g - inverse(jac) * f
    result = g1.mapIt(it[0])
    inc iter
    if iter > MaxIter: break
    var exit = true
    for idx, val in result:
      if abs(val - gu[idx]) > Tol:
        exit = false
        break
    if exit: break


when isMainModule:

  #[ Solve the two non-linear equations:
     y = -x^2 + x + 0.5
     y + 5xy = x^2
     given initial guesses of x = y = 1.2

     Example taken from:
     http://www.fixoncloud.com/Home/LoginValidate/OneProblemComplete_Detailed.php?problemid=286

     Expected results: x = 1.23332, y = 0.2122
  ]#

  let
    f1: Func = (x: Vector) => -x[0] * x[0] + x[0] + 0.5 - x[1]
    f2: Func = (x: Vector) => x[1] + 5 * x[0] * x[1] - x[0] * x[0]
    funcs1: Funcs = @[f1, f2]
    jacobian1: Jacobian = @[@[Func((x: Vector) => - 2 * x[0] + 1),
                              Func((x: Vector) => -1.0)],
                            @[Func((x: Vector) => 5 * x[1] - 2 * x[0]),
                              Func((x: Vector) => 1 + 5 * x[0])]]

    guesses1 = @[1.2, 1.2]
    sol1 = solve(funcs1, jacobian1, guesses1)
  echo &"Approximate solutions are x = {sol1[0]:.7f},  y = {sol1[1]:.7f}"

  #[ Solve the three non-linear equations:
     9x^2 + 36y^2 + 4z^2 - 36 = 0
     x^2 - 2y^2 - 20z = 0
     x^2 - y^2 + z^2 = 0
     given initial guesses of x = y = 1.0 and z = 0.0

     Example taken from:
     http://mathfaculty.fullerton.edu/mathews/n2003/FixPointNewtonMod.html (exercise 5)

     Expected results: x = 0.893628, y = 0.894527, z = -0.0400893
  ]#

  echo()
  let
    f3: Func = (x: Vector) => 9 * x[0] * x[0] + 36 * x[1] * x[1] + 4 * x[2] * x[2] - 36
    f4: Func = (x: Vector) => x[0] * x[0] - 2 * x[1] * x[1] - 20 * x[2]
    f5: Func = (x: Vector) => x[0] * x[0] - x[1] * x[1] + x[2] * x[2]
    funcs2: Funcs = @[f3, f4, f5]
    jacobian2: Jacobian = @[@[Func((x: Vector) => 18 * x[0]),
                              Func((x: Vector) => 72 * x[1]),
                              Func((x: Vector) =>  8 * x[2])],
                            @[Func((x: Vector) =>  2 * x[0]),
                              Func((x: Vector) => -4 * x[1]),
                              Func((x: Vector) => -20.0)],
                            @[Func((x: Vector) =>  2 * x[0]),
                              Func((x: Vector) => -2 * x[1]),
                              Func((x: Vector) =>  2 * x[2])]]
    guesses2 = @[1.0, 1.0, 0.0]
    sol2 = solve(funcs2, jacobian2, guesses2)
  echo &"Approximate solutions are x = {sol2[0]:.7f},  y = {sol2[1]:.7f}, z = {sol2[2]:.7f}"</syntaxhighlight>

{{out}}
<pre>Approximate solutions are x = 1.2333178,  y = 0.2122450

Approximate solutions are x = 0.8936282,  y = 0.8945270, z = -0.0400893</pre>

=={{header|Phix}}==
{{trans|Go}}
Uses code from [[Reduced_row_echelon_form#Phix]],
[[Gauss-Jordan_matrix_inversion#Phix]],
[[Matrix_transposition#Phix]], and 
[[Matrix_multiplication#Phix]]<br>
See std distro for a complete runnable version.
<!--<syntaxhighlight lang="phix">(phixonline)-->
 <span style="color: #000080;font-style:italic;">-- demo\rosetta\Multidimensional_Newton-Raphson_method.exw</span>
 <span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">solve</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">fs</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">jacob</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">guesses</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">integer</span> <span style="color: #000000;">size</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fs</span><span style="color: #0000FF;">),</span>
             <span style="color: #000000;">maxIter</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">12</span><span style="color: #0000FF;">,</span>
             <span style="color: #000000;">iter</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
     <span style="color: #004080;">sequence</span> <span style="color: #000000;">gu1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">g</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">g1</span><span style="color: #0000FF;">,</span>
              <span style="color: #000000;">gu2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">guesses</span><span style="color: #0000FF;">,</span>
              <span style="color: #000000;">jac</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">size</span><span style="color: #0000FF;">),</span><span style="color: #000000;">size</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">atom</span> <span style="color: #000000;">tol</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1e-8</span>
     <span style="color: #008080;">while</span> <span style="color: #004600;">true</span> <span style="color: #008080;">do</span>
         <span style="color: #000000;">gu1</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">gu2</span>
         <span style="color: #000000;">g</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">matrix_transpose</span><span style="color: #0000FF;">({</span><span style="color: #000000;">gu1</span><span style="color: #0000FF;">})</span>
         <span style="color: #000000;">t</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">size</span><span style="color: #0000FF;">)</span>
         <span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">size</span> <span style="color: #008080;">do</span>
             <span style="color: #000000;">t</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">call_func</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fs</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],{</span><span style="color: #000000;">gu1</span><span style="color: #0000FF;">})</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
         <span style="color: #000000;">f</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">matrix_transpose</span><span style="color: #0000FF;">({</span><span style="color: #000000;">t</span><span style="color: #0000FF;">})</span>
         <span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">size</span> <span style="color: #008080;">do</span>
             <span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">size</span> <span style="color: #008080;">do</span>
                 <span style="color: #000000;">jac</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">call_func</span><span style="color: #0000FF;">(</span><span style="color: #000000;">jacob</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">],{</span><span style="color: #000000;">gu1</span><span style="color: #0000FF;">})</span>
             <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
         <span style="color: #000000;">g1</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g</span><span style="color: #0000FF;">,</span><span style="color: #000000;">matrix_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inverse</span><span style="color: #0000FF;">(</span><span style="color: #000000;">jac</span><span style="color: #0000FF;">),</span><span style="color: #000000;">f</span><span style="color: #0000FF;">))</span>
         <span style="color: #000000;">gu2</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">vslice</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
         <span style="color: #000000;">iter</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
         <span style="color: #004080;">bool</span> <span style="color: #000000;">any</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">find</span><span style="color: #0000FF;">(</span><span style="color: #004600;">true</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sq_gt</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gu2</span><span style="color: #0000FF;">),</span><span style="color: #000000;">gu1</span><span style="color: #0000FF;">),</span><span style="color: #000000;">tol</span><span style="color: #0000FF;">))!=</span><span style="color: #000000;">0</span>
         <span style="color: #008080;">if</span> <span style="color: #008080;">not</span> <span style="color: #000000;">any</span> <span style="color: #008080;">or</span> <span style="color: #000000;">iter</span> <span style="color: #0000FF;">>=</span> <span style="color: #000000;">maxIter</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
     <span style="color: #008080;">return</span> <span style="color: #000000;">gu2</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 
 <span style="color: #008080;">function</span> <span style="color: #000000;">f1</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">-</span><span style="color: #000000;">y</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">f2</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">+</span><span style="color: #000000;">5</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">y</span><span style="color: #0000FF;">-</span><span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">f3</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">36</span><span style="color: #0000FF;">*</span><span style="color: #000000;">y</span><span style="color: #0000FF;">*</span><span style="color: #000000;">y</span><span style="color: #0000FF;">+</span><span style="color: #000000;">4</span><span style="color: #0000FF;">*</span><span style="color: #000000;">z</span><span style="color: #0000FF;">*</span><span style="color: #000000;">z</span><span style="color: #0000FF;">-</span><span style="color: #000000;">36</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">f4</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">y</span><span style="color: #0000FF;">*</span><span style="color: #000000;">y</span><span style="color: #0000FF;">-</span><span style="color: #000000;">20</span><span style="color: #0000FF;">*</span><span style="color: #000000;">z</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">f5</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">y</span><span style="color: #0000FF;">*</span><span style="color: #000000;">y</span><span style="color: #0000FF;">+</span><span style="color: #000000;">z</span><span style="color: #0000FF;">*</span><span style="color: #000000;">z</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 
 <span style="color: #008080;">function</span> <span style="color: #000000;">j1</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">j2</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000080;font-style:italic;">/*v*/</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">j3</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">*</span><span style="color: #000000;">y</span><span style="color: #0000FF;">-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">j4</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">5</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">j11</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #000000;">18</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">j12</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{?,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #000000;">72</span><span style="color: #0000FF;">*</span><span style="color: #000000;">y</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">j13</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{?,?,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #000000;">8</span><span style="color: #0000FF;">*</span><span style="color: #000000;">z</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">j21</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">j22</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{?,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">4</span><span style="color: #0000FF;">*</span><span style="color: #000000;">y</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">j23</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000080;font-style:italic;">/*v*/</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">20</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">j31</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">j32</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{?,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">y</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">j33</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{?,?,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span> <span style="color: #008080;">return</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">z</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 
 <span style="color: #008080;">procedure</span> <span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
 <span style="color: #004080;">sequence</span> <span style="color: #000000;">fs</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">jacob</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">guesses</span>
     <span style="color: #000080;font-style:italic;">/*
        solve the two non-linear equations:
        y = -x^2 + x + 0.5
        y + 5xy = x^2
        given initial guesses of x = y = 1.2
  
        Example taken from:
        http://www.fixoncloud.com/Home/LoginValidate/OneProblemComplete_Detailed.php?problemid=286
  
        Expected results: x = 1.23332, y = 0.2122
     */</span>
     <span style="color: #000000;">fs</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">}</span>
     <span style="color: #000000;">jacob</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">j1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j2</span><span style="color: #0000FF;">},</span>
              <span style="color: #0000FF;">{</span><span style="color: #000000;">j3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j4</span><span style="color: #0000FF;">}}</span>
     <span style="color: #000000;">guesses</span> <span style="color: #0000FF;">:=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1.2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.2</span><span style="color: #0000FF;">}</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Approximate solutions are x = %.7f,  y = %.7f\n\n"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">solve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fs</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">jacob</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">guesses</span><span style="color: #0000FF;">))</span>
  
     <span style="color: #000080;font-style:italic;">/*
        solve the three non-linear equations:
        9x^2 + 36y^2 + 4z^2 - 36 = 0
        x^2 - 2y^2 - 20z = 0
        x^2 - y^2 + z^2 = 0
        given initial guesses of x = y = 1.0 and z = 0.0
  
        Example taken from:
        http://mathfaculty.fullerton.edu/mathews/n2003/FixPointNewtonMod.html (exercise 5)
  
        Expected results: x = 0.893628, y = 0.894527, z = -0.0400893
     */</span>
  
     <span style="color: #000000;">fs</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">f3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">f4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">f5</span><span style="color: #0000FF;">}</span>
     <span style="color: #000000;">jacob</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">j11</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j12</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j13</span><span style="color: #0000FF;">},</span>
              <span style="color: #0000FF;">{</span><span style="color: #000000;">j21</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j22</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j23</span><span style="color: #0000FF;">},</span>
              <span style="color: #0000FF;">{</span><span style="color: #000000;">j31</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j32</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j33</span><span style="color: #0000FF;">}}</span>
     <span style="color: #000000;">guesses</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Approximate solutions are x = %.7f,  y = %.7f,  z = %.7f\n"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">solve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fs</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">jacob</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">guesses</span><span style="color: #0000FF;">))</span>
 
 <span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
 <span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Approximate solutions are x = 1.2333178,  y = 0.2122450

Approximate solutions are x = 0.8936282,  y = 0.8945270,  z = -0.04008929
</pre>

=={{header|Raku}}==
(formerly Perl 6)
<syntaxhighlight lang="raku" line># Reference:
# https://github.com/pierre-vigier/Perl6-Math-Matrix
# Mastering Algorithms with Perl
# By Jarkko Hietaniemi, John Macdonald, Jon Orwant
# Publisher: O'Reilly Media, ISBN-10: 1565923987
# https://resources.oreilly.com/examples/9781565923980/blob/master/ch16/solve

use v6;

sub solve_funcs ($funcs, @guesses, $iterations, $epsilon) {
   my ($error_value, @values, @delta, @jacobian); my \ε = $epsilon;
   for 1 .. $iterations {
      for ^+$funcs { @values[$^i] = $funcs[$^i](|@guesses); }
      $error_value = 0;
      for ^+$funcs { $error_value += @values[$^i].abs }
      return @guesses if $error_value ≤ ε;
      for ^+$funcs { @delta[$^i] = -@values[$^i] }
      @jacobian = jacobian $funcs, @guesses, ε;
      @delta = solve_matrix @jacobian, @delta;
      loop (my $j = 0, $error_value = 0; $j < +$funcs; $j++) {
         $error_value += @delta[$j].abs ;
         @guesses[$j] += @delta[$j];
      }
      return @guesses if $error_value ≤ ε;
   }
   return @guesses;
}

sub jacobian ($funcs is copy, @points is copy, $epsilon is copy) {
   my ($Δ, @P, @M, @plusΔ, @minusΔ);
   my Array @jacobian; my \ε = $epsilon;
   for ^+@points -> $i {
      @plusΔ = @minusΔ = @points;
      $Δ = (ε * @points[$i].abs) || ε;
      @plusΔ[$i] = @points[$i] + $Δ;
      @minusΔ[$i] = @points[$i] - $Δ;
      for ^+$funcs { @P[$^k] = $funcs[$^k](|@plusΔ); }
      for ^+$funcs { @M[$^k] = $funcs[$^k](|@minusΔ); }
      for ^+$funcs -> $j {
         @jacobian[$j][$i] = (@P[$j] - @M[$j]) / (2 * $Δ);
      }
   }
   return @jacobian;
}

sub solve_matrix (@matrix_array is copy, @delta is copy) {
   # https://github.com/pierre-vigier/Perl6-Math-Matrix/issues/56
   { use Math::Matrix;
      my $matrix = Math::Matrix.new(@matrix_array);
      my $vector = Math::Matrix.new(@delta.map({.list}));
      die "Matrix is not invertible" unless $matrix.is-invertible;
      my @result = ( $matrix.inverted dot $vector ).transposed;
      return @result.split(" ");
   }
}

my $funcs = [
   { 9*$^x² + 36*$^y² + 4*$^z² - 36 }
   { $^x² - 2*$^y² - 20*$^z }
   { $^x² - $^y² + $^z² }
];

my @guesses = (1,1,0);

my @solution = solve_funcs $funcs, @guesses, 20, 1e-8;

say "Solution: ", @solution;
</syntaxhighlight>
{{out}}
<pre>
Solution: [0.8936282344764825 0.8945270103905782 -0.04008928615915281]
</pre>

=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-matrix}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "./matrix" for Matrix
import "./fmt" for Fmt

var solve = Fn.new { |funcs, jacobian, guesses|
    var size = funcs.count
    var gu1 = []
    var gu2 = guesses.toList
    var jac = Matrix.new(size, size)
    var tol = 1e-8
    var maxIter = 12
    var iter = 0
    while (true) {
        gu1 = gu2
        var g = Matrix.new([gu1]).transpose
        var v = List.filled(size, 0)
        for (i in 0...size) v[i] = funcs[i].call(gu1)
        var f = Matrix.new([v]).transpose
        for (i in 0...size) {
            for (j in 0...size) jac[i, j] = jacobian[i][j].call(gu1)
        }
        var g1 = g - jac.inverse * f
        gu2 = List.filled(size, 0)
        for (i in 0...size) gu2[i] = g1[i][0]
        iter = iter + 1
        if (iter == maxIter) break
        if ((0...size).all { |i| (gu2[i] - gu1[i]).abs <= tol }) break
    }
    return gu2
}

/* solve the two non-linear equations:
   y = -x^2 + x + 0.5
   y + 5xy = x^2
   given initial guesses of x = y = 1.2

   Example taken from:
   http://www.fixoncloud.com/Home/LoginValidate/OneProblemComplete_Detailed.php?problemid=286

   Expected results: x = 1.23332, y = 0.2122
*/
var f1 = Fn.new { |x| -x[0] * x[0] + x[0] + 0.5 - x[1] }
var f2 = Fn.new { |x|  x[1] + 5 * x[0] * x[1] - x[0] * x[0] }
var funcs = [f1, f2]
var jacobian = [
    [ Fn.new { |x| - 2 * x[0] + 1 }, Fn.new { |x| -1 } ],
    [ Fn.new { |x| 5 * x[1] - 2 * x[0] }, Fn.new { |x| 1 + 5 * x[0] } ]
]
var guesses = [1.2, 1.2]
var sols = solve.call(funcs, jacobian, guesses)
Fmt.print("Approximate solutions are x = $.7f, y = $.7f", sols[0], sols[1])

/* solve the three non-linear equations:
   9x^2 + 36y^2 + 4z^2 - 36 = 0
   x^2 - 2y^2 - 20z = 0
   x^2 - y^2 + z^2 = 0
   given initial guesses of x = y = 1.0 and z = 0.0

   Example taken from:
   http://mathfaculty.fullerton.edu/mathews/n2003/FixPointNewtonMod.html (exercise 5)

   Expected results: x = 0.893628, y = 0.894527, z = -0.0400893
*/
System.print()
var f3 = Fn.new { |x| 9 * x[0] * x[0] + 36 * x[1] *x[1] + 4 * x[2] * x[2] - 36 }
var f4 = Fn.new { |x| x[0] * x[0] - 2 * x[1] * x[1] - 20 * x[2] }
var f5 = Fn.new { |x| x[0] * x[0] - x[1] * x[1] + x[2] * x[2] }
funcs = [f3, f4, f5]
jacobian = [
    [ Fn.new { |x| 18 * x[0] }, Fn.new { |x| 72 * x[1] }, Fn.new { |x| 8 * x[2] } ],
    [ Fn.new { |x|  2 * x[0] }, Fn.new { |x| -4 * x[1] }, Fn.new { |x| -20 } ],
    [ Fn.new { |x|  2 * x[0] }, Fn.new { |x| -2 * x[1] }, Fn.new { |x| 2 * x[2] } ]
]
guesses = [1, 1, 0]
sols = solve.call(funcs, jacobian, guesses)
Fmt.print("Approximate solutions are x = $.7f, y = $.7f, z = $.7f", sols[0], sols[1], sols[2])</syntaxhighlight>

{{out}}
<pre>
Approximate solutions are x = 1.2333178, y = 0.2122450

Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.0400893
</pre>

=={{header|zkl}}==
This doesn't use Newton-Raphson (with derivatives) but a hybrid algorithm.
<syntaxhighlight lang="zkl">var [const] GSL=Import.lib("zklGSL");    // libGSL (GNU Scientific Library)

   // two functions of two variables: f(x,y)=0
fs:=T(fcn(x,y){ x.atan() - y*y*y }, fcn(x,y){ 4.0*x*x + 9*y*y - 36 });
v=GSL.VectorFromData(2.75, 1.25);	// an initial guess at the solution
GSL.multiroot_fsolver(fs,v);
v.format(11,8).println();		// answer overwrites initial guess

fs.run(True,v.toList().xplode()).println();	// deltas from zero</syntaxhighlight>
{{out}}
<pre>
 2.59807621, 1.06365371
L(2.13651e-09,2.94321e-10)
</pre>
A condensed solver (for a different set of functions):
<syntaxhighlight lang="zkl">v:=GSL.VectorFromData(-10.0, -15.0);
GSL.multiroot_fsolver(T( fcn(x,y){ 1.0 - x }, fcn(x,y){ 10.0*(y - x*x) }),v)
.format().println();	// --> (1,1)</syntaxhighlight>
{{out}}
<pre>
1.00,1.00
</pre>
Another example:
<syntaxhighlight lang="zkl">v:=GSL.VectorFromData(1.0, 1.0, 0.0);	// initial guess
fxyzs:=T(
   fcn(x,y,z){ x*x*9 + y*y*36 + z*z*4 - 36 }, // 9x^2 + 36y^2 + 4z^2 - 36 = 0
   fcn(x,y,z){ x*x - y*y*2 - z*20 },	      // x^2 - 2y^2 - 20z = 0
   fcn(x,y,z){ x*x - y*y + z*z });	      // x^2 - y^2 + z^2 = 0
(v=GSL.multiroot_fsolver(fxyzs,v)).format(12,8).println();

fxyzs.run(True,v.toList().xplode()).println();	// deltas from zero</syntaxhighlight>
{{out}}
<pre>
  0.89362824,  0.89452701, -0.04008929
L(6.00672e-08,1.0472e-08,9.84017e-09)
</pre>