Capture-avoiding substitution function -- Lambda calculus

I am trying to write a function that performs capture-avoiding substitution in Lambda calculus. The code compiles but does not spit out the correct answer. I've written what I expect the code to do, is my comprehension correct?

For example, I should get the following output for this input (numeral 0 is the Church numeral 0)

*Main> substitute "b" (numeral 0) example    -- \a. \x. ((\y. a) x) b
\c. \a. (\a. c) a (\f. \x. x)

-- The incorrect result I actually got
\c. \c. (\f. \x. x) (x (\b. a))

NB \y is renamed to \a due to the substitution (\y.a)[N/b] (I think I have this covered in the code I have written, but please let me know if I am wrong.)

import Data.Char
import Data.List

type Var = String

data Term =
    Variable Var
  | Lambda   Var  Term
  | Apply    Term Term
  --  deriving Show

instance Show Term where
  show = pretty

example :: Term        -- \a. \x. ((\y. a) x) b
example = Lambda "a"
            (Lambda "x" (Apply (Apply (Lambda "y" (Variable "a")) 
                                      (Variable "x")) 
                               (Variable "b")))

pretty :: Term -> String
pretty = f 0
    where
      f i (Variable x) = x
      f i (Lambda x m) = if i /= 0 then "(" ++ s ++ ")" else s 
                         where s = "\\" ++ x ++ ". " ++ f 0 m 
      f i (Apply  n m) = if i == 2 then "(" ++ s ++ ")" else s 
                         where s = f 1 n ++ " " ++ f 2 m

substitute :: Var -> Term -> Term -> Term

substitute x n (Variable y)  
    --if y = x, then leave n alone   
    | y == x    = n
    -- otherwise change to y  
    | otherwise = Variable y

substitute x n (Lambda y m)
    --(\y.M)[N/x] = \y.M if y = x 
    | y == x    = Lambda y m
    --otherwise \z.(M[z/y][N/x]), where `z` is a fresh variable name 
    --generated by the `fresh` function, `z` must not be used in M or N, 
    --and `z` cannot be equal `x`. The `used` function checks if a 
    --variable name has been used in `Lambda y m`   
    | otherwise = Lambda newZ newM
                  where newZ = fresh(used(Lambda y m))
                        newM = substitute x n m          

substitute x n (Apply  m2 m1) = Apply newM2 newM1
    where newM1 = substitute x n m2
          newM2 = substitute x n m1

used :: Term -> [Var]
used (Variable n) = [n]
used (Lambda n t) = merge [n] (used t)
used (Apply t1 t2) = merge (used t1) (used t2)

variables :: [Var]
variables =  [l:[] | l <- ['a'..'z']] ++ 
             [l:show x | x <- [1..], l <- ['a'..'z']]

filterFreshVariables :: [Var] -> [Var] -> [Var]
filterFreshVariables lst = filter ( `notElem` lst)

fresh :: [Var] -> Var
fresh lst = head (filterFreshVariables lst variables)

recursiveNumeral :: Int -> Term
recursiveNumeral i
  | i == 0 = Variable "x"
  | i > 0 = Apply(Variable "f")(recursiveNumeral(i-1))

numeral :: Int -> Term
numeral i = Lambda "f" (Lambda "x" (recursiveNumeral i))

merge :: Ord a => [a] -> [a] -> [a]
merge (x : xs) (y : ys)
  | x < y = x : merge xs (y : ys)
  | otherwise = y : merge (x : xs) ys
merge xs [] = xs
merge [] ys = ys

Solution

This part in substitute x n (Lambda y m) is not correct:

the comment says "z must not be used in M or N", but there is nothing preventing that. newZ could be a variable in n, which leads to a problematic capture
the substitution z/y has not been done

    | otherwise = Lambda newZ newM
                  where newZ = fresh(used(Lambda y m))
                        newM = substitute x n m

Fix:

"z must not be used in M or N":

newZ = fresh(used m `merge` used n)

"M[z/y][N/x]":

newM = substitute x n (substitute y (Variable newZ) m)

Put together:

    | otherwise = Lambda newZ newM
    where
      newZ = fresh(used m `merge` used n)
      newM = substitute x n (substitute y (Variable newZ) m)

Note that refreshing all bindings as done above makes it difficult to understand the result and to debug substitution. Actually y only needs to be refreshed if y is in n. Otherwise you can keep y, adding this clause:

    | y `notElem` used n = Lambda y (substitute x n m)

Another idea would be to modify fresh to pick a name similar to the old one, e.g., by appending numbers until one doesn't clash.

There is still a bug I missed: newZ should also not be equal to x (the variable originally being substituted).

-- substitute [a -> \f. \x. x] in (\g. g), should be (\g. g)
ghci> substitute "a" (numeral 0) (Lambda "g" (Variable "g"))
\a. \g. \x. x

Two ways to address this:

add x to the set of variables to exclude newZ from:

newZ = fresh ([x] `merge` used m `merge` used n)

if you think about it, this bug only manifests itself when x is not in m, in which case there is nothing to substitute, so another way is to add one more branch skipping the work:
```
| x `notElem` used m = Lambda y m
```

Put together:

substitute x n (Lambda y m)
    --(\y.M)[N/x] = \y.M if y = x 
    | y == x    = Lambda y m
    | x `notElem` used m = Lambda y m
    | y `notElem` used n = Lambda y (substitute x n m)
    | otherwise = Lambda newZ newM
                  where newZ = fresh(used m `merge` used n)
                        newM = substitute x n (substitute y (Variable newZ) m)

Output

ghci> example
\a. \x. (\y. a) x b
ghci> numeral 0
\f. \x. x
ghci> substitute "b" (numeral 0) example
\a. \c. (\y. a) c (\f. \x. x)

Note: I haven't tried to prove this code correct (exercise for the reader: define "correct"), there may still be bugs I missed. There must be some course about lambda calculus that has all the details and pitfalls but I haven't bothered to look.