function y=myran(x); 
  a=1664525;
  b=1013904223;
  m=2^32;
  y=mod(a*x+b,m);
  return; 
end

x(1)=123;
for i=2:10; 
	x(i)=myran(x(i-1));
	disp(x(i)/2^32);
end

disp('Zera funkcji, metodami dwoma')
clear all % na wszelki wypadek, kasujemy wszystkie zmienne w octave
 
% definiujemy funkcje, ktore oblicza kolejne iterajec metod: Newtona i sieczncyh
function x=newton(x0,N,f,df)
	disp('Metoda Newtona')
	x=zeros(1,N);
	x(1)=x0;
	for i=2:N
		x(i)=x(i-1)-f(x(i-1))/df(x(i-1));
	endfor
endfunction
function x=sieczna(x0,x1,N,f)
	disp('Metoda siecznych')
	x=zeros(1,N);
	x(1)=x0;
	x(2)=x1;
	for i=3:N
		s = (f(x(i-1)) - f(x(i-2)) ) /(x(i-1)-x(i-2));
		x(i)=x(i-1)-f(x(i-1))/s;
	endfor
endfunction
 
% to jest funkcja pomocnicza, wykonujaca wykres
function plotall(newt,siecz,f,a,b)
	subplot(2,1,1)
	plot(newt,"r*-")
	hold on
	plot(siecz,"b*-")
	hold off
	subplot(2,1,2)
	x=linspace(a,b,101);
	plot(x,f(x))
	hold on
	plot(newt,f(newt),"r*")
	plot(siecz,f(siecz),"g*")
	hold off
endfunction
 
% teraz mozemy zdefiniowac funkcje matemetyczne, np f,g
% i oblicznyc ich zera  
f=@(x)  sign(x).*sqrt(abs(x));
df=@(x) 1/(2*sqrt(abs(x)));
N=10 
newt=newton(1.1,N,f,df);
siecz=sieczna(1.5,0.4,N,f,df);
figure(1) % figure pozwoli otworzyc kilka okien z wykresami
plotall(newt,siecz,f,-2.1,2)
 
g=@(x) sin(x)./x;
dg=@(x) cos(x)./x-sin(x)./(x.^2);
N=10
newt=newton(1.1,N,g,dg);
siecz=sieczna(1.5,0.4,N,g);
figure(2)
plotall(newt,siecz,g,0.01,12)
 
h=@(x)  atan(x);
dh=@(x) 1.0./(1+x.^2);
N=5
newt=newton(1.1,N,h,dh);
siecz=sieczna(1.5,0.4,N,h);
figure(3)
plotall(newt,siecz,h,-3.001,4)