function newtonraphson

stupen=2;
h=0.01;

function f1=f1(x,y)
    f1=x^2 + 4*x - y^2 - 2*y - 1;
end

function df1dx=df1dx(x,y,h)
    df1dx= 1/(2*h) * ( f1(x+h,y) - f1(x-h,y) );
end

function df1dy=df1dy(x,y,h)
    df1dy=1/(2*h) * ( f1(x,y+h) - f1(x,y-h) );
end

function f2=f2(x,y)
    f2=x^2 + 5*y - 4;
end

function df2dx=df2dx(x,y,h)
    df2dx = 1/(2*h) * ( f2(x+h,y) - f2(x-h,y) );
end

function df2dy=df2dy(x,y,h)
    df2dy = 1/(2*h) * ( f2(x,y+h) - f2(x,y-h) );
end

n=20;
x0=1;
y0=1;

a=zeros(stupen);
b=zeros(stupen,1);

for i=1:n

	disp(['Krok ' num2str(i) ' (x0,y0)=(' num2str(x0) ',' num2str(y0) ...
	')   (f1,f2)=(' num2str(f1(x0,y0)) ',' ...
	num2str(f2(x0,y0)) ')' ] );

	a(1,1)=df1dx(x0,y0,h);
	a(1,2)=df1dy(x0,y0,h);
	a(2,1)=df2dx(x0,y0,h);
	a(2,2)=df2dy(x0,y0,h);

	b(1) = -f1(x0,y0);
	b(2) = -f2(x0,y0);

	dx=a\b;

	x0=x0+dx(1);
	y0=y0+dx(2);

end

end