Метод Ньютона безусловно хорош, в случае СНАУ думаю это лучший выбор. Правда у него часто бывает очень узкая область сходимости. Но это решаемо, вводится настраиваемый параметр сходимости
![$\theta[k+1]=\theta[k]-\mu J^{-1}\Delta$ $\theta[k+1]=\theta[k]-\mu J^{-1}\Delta$](https://dxdy-01.korotkov.co.uk/f/c/f/8/cf8c4828a047dab9fe41da8e20ad2bd182.png)
, где
![$J, \Delta$ $J, \Delta$](https://dxdy-03.korotkov.co.uk/f/e/c/c/ecc23aea2d53dca7b57cbb8f618f7c9782.png)
- матрица Якоби и невязка.
При
![$\mu=1$ $\mu=1$](https://dxdy-03.korotkov.co.uk/f/2/c/c/2cc5b0184fa1062f086c02d06101eabd82.png)
это классический метод Ньютона, при уменьшении этого параметра область сходимости расширяется, правда сходимость становится линейной. Поэтому по мери приближения к решению
![$\mu$ $\mu$](https://dxdy-01.korotkov.co.uk/f/0/7/6/07617f9d8fe48b4a7b3f523d6730eef082.png)
желательно постепенно уменьшать до нуля, разумеется при сохранении сходимости.
Если пишите на С, то думаю связываться с каким либо специализированным ПО не стоит, потратите больше времени на его изучение, да и результаты потом не сможете использовать как вам захочется (конечно если это не одноразовое вычисление).