diff --git a/hafar.cabal b/hafar.cabal
--- a/hafar.cabal
+++ b/hafar.cabal
@@ -4,10 +4,10 @@
 --
 -- see: https://github.com/sol/hpack
 --
--- hash: 11f7d969e8f00b5e392eda2b0897debe20363e72178cbd6a8a87a7ab29e2b819
+-- hash: 6bf5e0c8d0a1cae72591bdd92749728d178e27c39856dbdc12dc983ea2f322f6
 
 name:           hafar
-version:        0.1.0.0
+version:        0.1.1.0
 synopsis:       Affine arithmetic library for Haskell
 description:    Hafar is an affine arithmetic library for Haskell. It is an efficient way to work with ranges of values or imprecise values.
 category:       Numeric
@@ -31,13 +31,14 @@
   exposed-modules:
       Numeric.AffineForm
       Numeric.AffineForm.ExplicitRounding
+      Numeric.AffineForm.Subdivision
   other-modules:
       Numeric.AffineForm.Utils
       Numeric.AffineForm.Internal
   hs-source-dirs:
       src
   build-depends:
-      base >=4.12 && <4.14
+      base >=4.12 && <4.13
     , intervals >=0.8 && <0.9
     , mtl >=2.2 && <2.3
   default-language: Haskell2010
@@ -50,13 +51,14 @@
       Numeric.AffineForm.Internal
       Numeric.AffineForm.Utils
       Numeric.AffineForm.ExplicitRounding
+      Numeric.AffineForm.Subdivision
   hs-source-dirs:
       test
       src
   ghc-options: -threaded -rtsopts -with-rtsopts=-N
   build-depends:
       QuickCheck >=2.13 && <2.14
-    , base >=4.12 && <4.14
+    , base >=4.12 && <4.13
     , hafar
     , intervals >=0.8 && <0.9
     , mtl >=2.2 && <2.3
diff --git a/src/Numeric/AffineForm.hs b/src/Numeric/AffineForm.hs
--- a/src/Numeric/AffineForm.hs
+++ b/src/Numeric/AffineForm.hs
@@ -6,12 +6,13 @@
                            midpoint,
                            inf, sup,
                            interval,
+                           afError,
                            member,
                            epscount_,
                            setMidpoint,
                            fix,
                            addError,
-                           (.+), (.*)
+                           (.+), (.*),
                           ) where
 
 import Numeric.AffineForm.Internal
diff --git a/src/Numeric/AffineForm/Internal.hs b/src/Numeric/AffineForm/Internal.hs
--- a/src/Numeric/AffineForm/Internal.hs
+++ b/src/Numeric/AffineForm/Internal.hs
@@ -15,6 +15,7 @@
 
 import Numeric.AffineForm.Utils
 import Numeric.AffineForm.ExplicitRounding
+import Numeric.AffineForm.Subdivision
 import qualified Numeric.Interval as IA
 import Numeric.Interval ((...))
 import Data.Fixed (mod')
@@ -69,32 +70,35 @@
   acosh = minrange acosh (\x -> 1/((sqrt (x-1))*(sqrt (x+1)))) Concave
   atanh = minrange atanh (\x -> 1/(1-x^2)) undefined
 
+instance (Ord a, Fractional a, ExplicitRounding a) => Subdivisible (AF s a) where
+  subdivide af n = mapToInterval af <$> subdivide (interval af) n
+  combine2 l r   = fromInterval $ combine2 (interval l) (interval r)
+
 type AFIndex = Int
 
 -- | AFM is a state monad that ensures that any new noise symbols have not been used by any previous affine form.
 -- All affine arithmetic calculations should be done inside the AFM monad. Affine forms do not make sense outside of their monad context.
-newtype AFMT t s m a = AFMT {runAFMT :: s -> m (a, s)}
-type AFM t a = AFMT t AFIndex Identity a
+newtype AFMT t m a = AFMT {runAFMT :: AFIndex -> m (a, AFIndex)}
+type AFM t a = AFMT t Identity a
 
-instance (Monad m) => Functor (AFMT t s m) where
+instance (Monad m) => Functor (AFMT t m) where
   fmap = liftM
 
-instance (Monad m) => Applicative (AFMT t s m) where
+instance (Monad m) => Applicative (AFMT t m) where
   pure = return
   (<*>) = ap
 
-instance (Monad m) => Monad (AFMT t s m) where
+instance (Monad m) => Monad (AFMT t m) where
   return a = AFMT $ \s -> return (a, s)
   (AFMT x) >>= f = AFMT $ \s -> do
     (v, s') <- x s
     (AFMT x') <- return $ f v
     x' s'
 
-instance (Monad m) => MonadState s (AFMT t s m) where
-  get   = AFMT $ \s -> return (s, s)
-  put s = AFMT $ \_ -> return ((), s)
+getIndex   = AFMT $ \s -> return (s, s)
+putIndex s = AFMT $ \_ -> return ((), s)
 
-instance MonadTrans (AFMT t s) where
+instance MonadTrans (AFMT t) where
   lift c = AFMT $ \s -> c >>= (\x -> return (x, s))
 
 -- | This gives an affine form with midpoint 0 and radius 1.
@@ -102,8 +106,8 @@
 -- It can be used to instantiate new affine forms.
 newEps :: Num a => AFM t (AF t a)
 newEps = do
-  idx <- get
-  put $ idx + 1
+  idx <- getIndex
+  putIndex $ idx + 1
   return $ AF 0 (replicate idx 0 ++ [1]) 0
 
 -- | Creates a new affine form that covers the interval.
@@ -119,12 +123,12 @@
 singleton x = AF x [] 0
 
 -- | Creates a new affine form that approximately represents some value.
--- This function adds a small error to account for the 'wobble' in the computer representation of the value.
+-- This function adds a small error to account for the `wobble` in the computer representation of the value.
 approxSingleton :: (ExplicitRounding a) => a -> AF s a
 approxSingleton x = AF x [] $ eps x
 
 -- | Evaluates the AFM monad. It is not possible to get an AF out of an AFM monad.
-evalAFM :: forall a b. (forall t. AFM t b) -> b
+evalAFM :: (forall t. AFM t b) -> b
 evalAFM (AFMT x) = fst . runIdentity $ x 0
 
 -- | Gives the radius of the affine form
@@ -158,6 +162,10 @@
 epscount_ :: AF s a -> Int
 epscount_ (AF _ xs _) = length xs
 
+-- | Returns the value of the error term
+afError :: AF s a -> a
+afError (AF _ _ e) = e
+
 -- Affine arithmetic operations
 
 -- | Sets the midpoint of the affine form
@@ -169,6 +177,15 @@
 addError (AF x xs xe) e
   | e >= 0 = AF x xs (xe + e)
   | otherwise = throw AddingNegativeError
+
+mapToInterval :: (Fractional a, ExplicitRounding a) => AF s a -> IA.Interval a -> AF s a
+mapToInterval af i = m .+ ((wi/waf) .* af)
+  where waf = radius af
+        wi  = IA.width i
+        m   = IA.midpoint i
+
+fromInterval :: (Fractional a) => IA.Interval a -> AF s a
+fromInterval i = AF (IA.midpoint i) [] ((IA.width i)/2)
 
 -- | Adds a scalar value to the affine form
 (.+) :: (Num a, ExplicitRounding a) => a -> AF s a -> AF s a
diff --git a/src/Numeric/AffineForm/Subdivision.hs b/src/Numeric/AffineForm/Subdivision.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AffineForm/Subdivision.hs
@@ -0,0 +1,125 @@
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE DeriveFunctor #-}
+
+-- | Subdivision contains an implementation of the branch-and-bound algorithm.
+-- This method divides interval-like values into smaller subdivisions and then applies a function to those values.
+-- This results in a more accurate output at the cost of having to repeat the calculations multiple times
+module Numeric.AffineForm.Subdivision
+  ( Subdivisible(..)
+  , SubdivisionEnvironment(..)
+  , defaultEnvironment
+  , branchAndBound
+  ) where
+
+import Control.Monad.Reader
+import Control.Monad.State
+
+import Data.List
+
+import Numeric.Interval as IA
+
+-- | The 'Subdivisible' class is used for datatypes that can be broken down into smaller pieces
+-- and combined together
+--
+-- The 'subdivide' function subdivides the value into `n` smaller values and 'combine' joins
+-- the subdivisions together into a bigger value
+--
+-- Subdividing and then combining a value is expected to give the initial value (sans rounding errors)
+class Subdivisible a where
+  subdivide    :: a -> Int -> [a]
+  combine2     :: a -> a -> a
+  combine      :: [a] -> a
+
+  combine2 l r = combine [l,r]
+  combine      = foldl1 combine2
+
+instance (Ord a, Fractional a) => Subdivisible (Interval a) where
+  subdivide i n = f . fromIntegral <$> [1..n]
+    where f x = (((x-1)*w/n')IA....(x*w/n')) + (IA.singleton $ inf i)
+          w   = width i
+          n'  = fromIntegral n
+  combine2 = hull
+
+instance (Subdivisible a) => Subdivisible [a] where
+  subdivide l n = sequenceA $ flip subdivide n <$> l
+  combine2 l r  = uncurry combine2 <$> zip l r
+
+-- | A data structure for configuring the 'branchAndBound' method
+--
+-- 'function' is the function that gets evaluated
+-- 'errorFun' is the error measuring function of the result
+-- 'maxError' specifies the maximum permitted error of an evaluation
+-- 'maxDepth' specifies how deep the subdivision can go
+-- 'subdivs'  specifies the number of subdivisions per value
+data SubdivisionEnvironment a b e = SubdivisionEnvironment
+  { function :: a -> b
+  , errorFun :: b -> e
+  , maxError :: e
+  , maxDepth :: Int
+  , subdivs  :: Int
+  }
+
+-- | This function creates a simple subdivision configuration.
+-- It requires the evaluator function and an error measuring function as its parameters.
+defaultEnvironment :: (Fractional e) => (a -> b) -> (b -> e) -> SubdivisionEnvironment a b e
+defaultEnvironment f g = SubdivisionEnvironment
+  { function = f
+  , errorFun = g
+  , maxError = 0.1
+  , maxDepth = 3
+  , subdivs  = 2
+  }
+
+type Subdivider a b e = Reader (SubdivisionEnvironment a b e) (SubdivisionTree a)
+
+data SubdivisionTree a
+  = Branch [SubdivisionTree a]
+  | Node a
+  deriving (Show, Functor)
+
+deepen :: (Subdivisible a, Ord e) => SubdivisionTree a -> Int -> Subdivider a b e
+deepen (Node x) depth =
+  do
+    env <- ask
+    let res = function env $ x
+        err = errorFun env $ res
+        n   = subdivs env
+    if err <= maxError env
+       -- Accurate answer found
+       then return $ Node x
+       -- Subdivide and deepen
+       else deepen (Branch $ Node <$> subdivide x n) $ depth
+deepen (Branch l) depth =
+  do
+    env <- ask
+    if depth >= maxDepth env
+       then return $ Branch l
+       else Branch <$> (sequence $ flip deepen (depth + 1) <$> l)
+
+collapse :: (Subdivisible a) => SubdivisionTree a -> a
+collapse (Node a)   = a
+collapse (Branch l) = combine $ collapse <$> l
+
+evalTree :: (Subdivisible b) => Subdivider a b e -> SubdivisionEnvironment a b e -> b
+evalTree div cfg = flip runReader cfg $
+  do
+    val <- div
+    env <- ask
+    return . collapse $ function env <$> val
+
+-- | This function iteratively subdivides a value to arrive at a more accurate result
+branchAndBound :: (Subdivisible a, Subdivisible b, Ord e) => a -> SubdivisionEnvironment a b e -> b
+branchAndBound x env = flip evalTree env $ deepen (Node x) 0
+
+test :: Interval Double
+test =
+  flip evalTree env $ deepen (Node $ [3 IA.... 4, 1 IA....4]) 0
+  where env = SubdivisionEnvironment
+          { function = (\[x, y] -> x * y - x)
+          , errorFun = width
+          , maxError = 0.1
+          , maxDepth = 5
+          , subdivs  = 2
+          }
+