Rozwiązywanie układu równań nieliniowych metodą itera

IDEA:
Dla uproszczenia rozważę układ dwóch równań. Niech będzie dany układ :
F1(x,y)=0 i F2(x,y)=0 (1)
Niech x0 i y0 są przybliżonymi pierwiastkami układu równań (1).
Przekształcając ten układ do postaci:
x=fi1(x,y) i y=fi2(x,y) można zbudować ciągi iteracyjne:

x_1=fi1(x_0,y_0) i y_1=fi2(x_0,y_0)
x_2=fi1(x_1,y_1) i y_2=fi2(x_1,y_1)
x_3=fi1(x_2,y_2) i y_3=fi2(x_2,y_2)
. .
. .
. .
x_n=fi1(x_(n-1),y_(n-1)) i y_4=fi2(x_(n-1),y_(n-1))
Uwaga!! Dla uproszczenia stosuje dwa oznaczenia indeksów dolnych (x1=x_1).

To już właściwie koniec, jeżeli proces iteracyjny jest zbieżny to przy dostatecznie dużym n różnice
|x_n-x_(n-1)|<eps oraz |y_n-y_(n-1)|<eps gdzie eps jest założoną dokładnością wyznaczania pierwiastków. Wtedy zachodzi, że:
F1(x_n,y_x) jest bliskie zeru i F2(x_n,y_n) jest bliskie zeru

Czyli chodzi o wyrażenie f(x) przez x=... i zastosowanie sie do powyższego algorytmu :)

Przykładowy program mojego autorstwa:
Uwaga! Crt2 to crt przystosowany do szybkich procesorów.
Program interakcja;
uses crt2;
const n=1300;
var x,y,x0,y0,dx,dy:double;
i:integer;
eps,epsp:double;

begin
clrscr;
write('Podaj x0, y0, eps: ');
read(x,y,eps);
i:=1;
x0:=sqrt((xy+5x-1)/2);
y0:=sqrt(x+3ln(x));
repeat
i:=i+1;
x:=sqrt((x0y0+5x0-1)/2);
y:=sqrt(x0+3ln(x0));
dx:=abs(x-x0);
dy:=abs(y-y0);
x0:=x;
y0:=y;
until (i>=1000)or((dx<=eps)and(dy<eps));
if i>=1000 then begin Write('Uklad nie ma rozwiazan');end;
writeln('x= ',x05,' ');
writeln('y= ',y05,' ');
writeln('F1(x)=',2x0x0-x0y0-5x0+13,'');
writeln('F2(x)=',x0+3ln(x0)-y0y03,'');
readkey;
end.

Oczywiście można tą metodę stosować to układów o wyszżym stopniu:)