One maybe not super-efficent but fun solution would be to abuse scipy.linalg.solve_banded
like so
import numpy as np
from scipy import linalg
N = 50
a = np.zeros((N,)) + np.array([[1, -1, -1]]).T
b = np.zeros((N,))
b[0] = 1
linalg.solve_banded((2, 0), a, b)
# array([ 1.00000000e+00, 1.00000000e+00, 2.00000000e+00,
# 3.00000000e+00, 5.00000000e+00, 8.00000000e+00,
# 1.30000000e+01, 2.10000000e+01, 3.40000000e+01,
# 5.50000000e+01, 8.90000000e+01, 1.44000000e+02,
# 2.33000000e+02, 3.77000000e+02, 6.10000000e+02,
# 9.87000000e+02, 1.59700000e+03, 2.58400000e+03,
# 4.18100000e+03, 6.76500000e+03, 1.09460000e+04,
# 1.77110000e+04, 2.86570000e+04, 4.63680000e+04,
# 7.50250000e+04, 1.21393000e+05, 1.96418000e+05,
# 3.17811000e+05, 5.14229000e+05, 8.32040000e+05,
# 1.34626900e+06, 2.17830900e+06, 3.52457800e+06,
# 5.70288700e+06, 9.22746500e+06, 1.49303520e+07,
# 2.41578170e+07, 3.90881690e+07, 6.32459860e+07,
# 1.02334155e+08, 1.65580141e+08, 2.67914296e+08,
# 4.33494437e+08, 7.01408733e+08, 1.13490317e+09,
# 1.83631190e+09, 2.97121507e+09, 4.80752698e+09,
# 7.77874205e+09, 1.25862690e+10])
How this works is that we write F_0, F_1, F_2... as one vector F and the identities -F_{i-1} -F_i + F{i+1} = 0 as a matrix which is nicely banded and then solve for F. Note that this can adapted to other similar recurrences.