To draw a fractal, the program applies Newton's method to complex numbers in a region of the complex plane. It then colors each point in the region to indicate the root of the equation to which Newton's method leads for that point.

For example, the equation Z^4 - 1 has four roots (values of Z that make the equation zero): 1, -1, i, and -i.

To make the calculations easier to understand, the program uses the following Complex class.

Option Explicit

Public Re As Double
Public Im As Double

Public Function Clone() As Complex
    Set Clone = NewComplex(Re, Im)
End Function

Public Sub Initialize(ByVal new_re As Double, ByVal new_im _
    As Double)
    Re = new_re
    Im = new_im
End Sub

Public Function MagnitudeSquared() As Double
    MagnitudeSquared = Re * Re + Im * Im
End Function

Public Function Magnitude() As Double
    Magnitude = Sqr(MagnitudeSquared())
End Function

Public Function Times(ByVal x As Complex) As Complex
    Set Times = NewComplex(Re * x.Re - Im * x.Im, Re * x.Im _
        + Im * x.Re)
End Function

Public Function TimesComponents(ByVal new_re As Double, _
     ByVal new_im As
Double) As Complex
    Set TimesComponents = NewComplex(Re * new_re - Im * _
         new_im, Re * new_im
+ Im * new_re)
End Function

Public Function Plus(ByVal x As Complex) As Complex
    Set Plus = NewComplex(Re + x.Re, Im + x.Im)
End Function

Public Function PlusComponents(ByVal new_re As Double, _
     ByVal new_im As
Double) As Complex
    Set PlusComponents = NewComplex(Re + new_re, Im + _
        new_im)
End Function

Public Function DividedBy(ByVal denominator As Complex) As _
    Complex
Dim conjugate As Complex
Dim numerator As Complex

    Set conjugate = denominator.Clone
    conjugate.Im = -conjugate.Im

    Set denominator = denominator.Times(conjugate)
    Set numerator = Me.Times(conjugate)
    numerator.Re = numerator.Re / denominator.Re
    numerator.Im = numerator.Im / denominator.Re

    Set DividedBy = numerator
End Function

Public Function IsCloseTo(ByVal x As Complex) As Boolean
    IsCloseTo = (Abs(Re - x.Re) + Abs(Im - x.Im) < 0.001)
End Function

Public Function ToString() As String
    If Im < 0 Then
        ToString = Format$(Re) & " - " & Format$(Abs(Im)) & _
            "i"
    Else
        ToString = Format$(Re) & " + " & Format$(Im) & "i"
    End If
End Function

Public Function NewComplex(ByVal new_re As Double, ByVal _
     new_im As Double)
As Complex
Dim result As Complex

    Set result = New Complex
    result.Initialize new_re, new_im
    Set NewComplex = result
End Function

' The function.
' F(x) = x^NUM_ROOTS - 1.
Private Function F(ByVal x As Complex) As Complex
Dim result As Complex
Dim i As Integer

    Set result = x.Clone()   ' x.

    ' x^NUM_ROOTS.
    For i = 1 To NUM_ROOTS - 1
        Set result = result.Times(x)
    Next i

    ' x^NUM_ROOTS - 1.
    Set F = result.PlusComponents(-1, 0)
End Function

' The function's derivative.
' dFdx(x) = NUM_ROOTS * x^(NUM_ROOTS - 1).
Private Function dFdx(ByVal x As Complex) As Complex
Dim result As Complex
Dim i As Integer

    Set result = x.Clone()   ' x.

    ' x^(NUM_ROOTS - 1).
    For i = 1 To NUM_ROOTS - 2
        Set result = result.Times(x)
    Next i

    ' NUM_ROOTS * x^(NUM_ROOTS - 1).
    Set dFdx = result.TimesComponents(NUM_ROOTS, 0)
End Function

' Draw the fractal.
Private Sub DrawFractal()
Const CUTOFF As Single = 0.0000001  ' Work until the
    ' epsilon squared < this.
Const MAX_ITER As Integer = 200     ' At most this many
    ' iterations.

Dim clr As Integer
Dim x As Complex
Dim x0 As Complex
Dim epsilon As Complex
Dim dx As Double
Dim dy As Double
Dim wid As Double
Dim hgt As Double
Dim iter As Integer
Dim r As Integer
Dim i As Integer
Dim j As Integer

    DrawMode = vbCopyPen
    Line (0, 0)-(ScaleWidth, ScaleHeight), BackColor, BF
    MousePointer = vbHourglass
    DoEvents

    ' Adjust the aspect ratio.
    AdjustAspect

    ' dx is the change in the real part
    ' (X value) for C. dy is the change in the
    ' imaginary part (Y value).
    wid = ScaleWidth
    hgt = ScaleHeight
    dx = (VisibleXmax - VisibleXmin) / (wid - 1)
    dy = (VisibleYmax - VisibleYmin) / (hgt - 1)

    ' Calculate the values.
    Set x0 = NewComplex(Xmin, 0)
    For i = 1 To wid
        x0.Im = Ymin
        For j = 1 To hgt
            Set x = x0
            iter = 0
            Do
                iter = iter + 1
                If iter > MAX_ITER Then Exit Do
                Set epsilon = _
                    F(x).DividedBy(dFdx(x)).TimesComponents(-1, _
                    0)
                Set x = x.Plus(epsilon)
            Loop While epsilon.MagnitudeSquared() > CUTOFF

            ' Set which root this is.
            If iter > MAX_ITER Then
                clr = 0 ' Black.
            Else
                ' Find the root.
                clr = 1 ' White.
                For r = 0 To UBound(m_Roots)
                    If x.IsCloseTo(m_Roots(r)) Then
                        clr = r + 2 ' +2 to skip black and
                            ' white.
                        Exit For
                    End If
                Next r
            End If

            ' Set the pixel's color.
            PSet (i, j), m_Colors(clr)

            ' Move to the next point.
            x0.Im = x0.Im + dy
        Next j
        x0.Re = x0.Re + dx
        DoEvents
    Next i

    Refresh
    Picture = Image

    MousePointer = vbCrosshair
    DrawMode = vbInvert
End Sub

Unfortunately this code is quite slow because of all of the heavy-lifting it does manipulating Complex objects. The VB .NET version is much faster.

For more information on Newton's method, see Eric W. Weisstein's article Newton's Method from MathWorld--A Wolfram Web Resource.