subroutine Muller_complex(f, h, z0, z1, z2, eps, iter, rad, kod)

complex(8) f
real(8) h

complex(8) z0, z1, z2, z3, rad
complex(8) a, b, c
complex(8) f10, f02, f21
complex(8) b1, b2
complex(8) delta

! Metoda Mulller
integer width

if(kod /=0) then
write(1, '(a, t15, a, t38, a, t61, a)') ' Iteratia', 'z0', 'z1', 'z2'
write(1, '(a)') ' * z3'
endif

lnit =iter

iter =0

100	iter =iter +1
c =f(z2)


f10 =(f(z0)-f(z1))/(z0-z1)
f02 =(f(z2)-f(z0))/(z2-z0)
f21 =(f(z1)-f(z2))/(z1-z2)

b =f21 +f02 -f10

a =(f10 -f21)/(z0-z2)

!*
!write(1,*) 'a, b, c:', a, b ,c

!delta =b**2-4*a*c

!write(1,*) 'delta:', delta
!*

delta =cdsqrt(b**2-4*a*c)

b1 =b +delta
b2 =b -delta

if(abs(b1) >abs(b2)) then
z3 =z2 -2*c/b1

else
z3 =z2 -2*c/b2
endif


if(kod /=0) then
iwiter =width(iter)

write(1, '(/, 1x, I<iwiter>, 3x, \)') iter
write(1,*) z0, z1, z2 
write(1,*) '*', z3
endif

if(abs(z3-z2) <= eps) then
rad =z3
kod =0
return

else
call select(z0, z1, z2, z3)
endif

if(iter <lnit) then
goto 100
else
rad =z3
kod =1
!iter =lnit
endif

end

