%!PS
%%BoundingBox: 10 40 600 350
/doos { newpath dup 3 2 roll dup 6 1 roll exch moveto
        3 2 roll dup 4 1 roll exch lineto
        dup 4 1 roll lineto  exch lineto  closepath stroke } def
%
% 10 40 600 350 doos
%
0 20 translate  150 150 scale
0.006 setlinewidth
%
% X-axis and Y-axis:
newpath 0.2 0.27 moveto 3.75 0 rlineto  0.27 0.2 moveto 0 1.95 rlineto  stroke
%
0.27 0.27 translate
%
% Start & increments:
/x 0 def   /dx { x /x x 0.01 add def } def
%
% Function  y = sqrt(w)  with Newton-Rhapson.
%
% Theory:  y(n+1) = y(n) - f(y(n))/f'(y(n))
%
% Here: y = sqrt(w)  or  y^2 = w . Take f(y) = y^2 - w , then  f'(y) = 2.y
% ==>   y(n+1) = y(n) - (y(n)^2 - w)/(2.y(n)) = y(n) - y(n)/2 + w/y(n)/2
%
% Conclusion:  y(n+1) = 1/2.( y(n) + w / y(n) )
%              --------------------------------
/wortel
{ /y 0.5 def  /w exch def  w 0 le { /y 0 def }
{ 7 { /y  y w y div add 0.5 mul  def } repeat } ifelse y } def
%
  /fun { dup 0 le { } { wortel } ifelse } def
%
  newpath dx dup fun moveto
  361 { dx dup fun lineto } repeat stroke
%
0.002 setlinewidth newpath
%
/Times-Roman findfont 0.10 scalefont setfont
%
/lijnen { x dup fun lineto  0 x sub 0 rlineto } def
/vervolgens { /x x fun def x 0 lineto lijnen } def
/x 3.5 def x 0 moveto lijnen 8 { vervolgens } repeat
/x 0.1 def x 0 moveto lijnen 8 { vervolgens } repeat
stroke
%
newpath
%
3.5 0.05 sub -0.08 moveto (3.5) show
0.1 0.05 sub -0.08 moveto (0.1) show
%
0.95 -0.08 moveto (1) show
-0.08 0.96 moveto (1) show
-0.08 -0.08 moveto (O) show
%
/Times-Roman findfont 0.13 scalefont setfont
2.5 -0.09 moveto (x) show
-0.09 1.6 moveto (y) show
2.7 dup fun 0.12 sub moveto (y\50x\51 = sqrt\50x\51) show
%
showpage