I've found a few online programs to generate first and follow sets from Backus-Naur form expressed as yacc (bison) grammars. They haven't worked for any but tiny examples. This program generates the first sets for bigger grammars.
Code:
#! /usr/bin/python3
# expand first.py > p.py && chmod +x p.py && ./p.py < latexp2clean.y

'''
    This program produces `first' sets from a yacc grammar definition.
    Prepare the input from a yacc grammar
    1) change empty productions to EPSILON
    2) remove c source from the rules section.
    3) rules start following the %% line as one would expect.

    bash use:

        $ chmod +x first.py
        $ ./first.py grammar.yacc

    Sample input (taken from my test suite, not included)
    
FIRST(S)={a,b,c}
FIRST(A)={a,EPSILON}
FIRST(B)={b,EPSILON}
FIRST(C)={c}
FIRST(D)={d,EPSILON}
FIRST(E)={e,EPSILON}
%%
S : A B C D E ; 
A : a | EPSILON ; 
B : b | EPSILON ; 
C : c ; 
D : d | EPSILON ; 
E : e | EPSILON ; 
%%
FOLLOW(S)={$}
FOLLOW(A)={b,c}
FOLLOW(B)={c}
FOLLOW(C)={d,e,$}
FOLLOW(D)={e,$}
FOLLOW(E)={$}

'''

import sys, collections, pprint, re, pdb

EPSILON = 'EPSILON'
epsilon = [EPSILON]
sepsilon = set(epsilon)

class Rule:

    def __init__(self, s, split = re.compile(' *[:|] *')):
        self.str = s
        fields = split.split(s)
        self.name = fields[0]
        self.productions = [production.split() for production in fields[1:]]

    def __getitem__(self, n):
        return self.productions[n]

    def __len__(self):
        return len(self.productions)

    def __str__(self):
        return self.str

class Graph(collections.OrderedDict):

    def __init__(self, rules):
        for rule in rules:
            self[rule.name] = rule.productions
        for rule in rules: # terminals
            for production in rule.productions:
                if production != epsilon:
                    for word in production:
                        if word not in self:
                            self[word] = None

#    @memoize...
    def first(self, key):
#        pdb.set_trace()
        return self._first(key, set(), 0)

    def _first(self, key, traced, depth):
        if key in traced:
            return set()
        traced.add(key)
        productions = self[key]
        if productions is None:  # terminal
            s = {key}
        else:
#            print(depth, key, productions)
            s = set()
            emptyok = False
            for production in productions:
                if production != epsilon:
                    for rule in production:
                        firsty = self._first(rule, traced, depth + 1)
#                        print(depth, key, productions, production, rule, firsty)
                        s = s.union(firsty)
                        if EPSILON not in firsty:
                            break
                    else:
                        emptyok = True
            s = s - sepsilon
            if (epsilon in productions) or emptyok:
                s.add(EPSILON)
        return s

def get_bnf(source):
    bnf = []
    accept = False
    for line in source:
        s = line.strip()
        accept ^= s == '%%'
        if accept and s:
           bnf.append(s)
    if '%%' in bnf[0]:
        del bnf[0]
#    pdb.set_trace()
    s = ''.join(bnf)
    s = (s[:-1] if s[-1] == ';' else s).strip()
    return list(map(Rule, s.split(';')))

def construct_graph(source):
    bnf = get_bnf(source)
    graph = Graph(bnf)
    return graph

def main(source = sys.stdin):
    graph = construct_graph(source)
    pprint.pprint(graph)
#    pdb.set_trace()
    for key in graph:
        first = graph.first(key)
        if first:
            print((key, first))

if __name__ == '__main__':
    if len(sys.argv) < 2:
        main()
    else:
        for filename in sys.argv[1:]:
            with open(filename, 'rt') as source:
                main(source)