Solving implicit equation with scipy.optimize.fsolve

I came to ask this question from the very similar one where I "learned" how to employ fsolve to solve for solving implicit equations.
My equation is defined by this function:

def fn4(Fs):
    return ( np.sum( (lamL * c + lamW * np.tan(np.radians(phi))) 
                 / (np.cos(lamA) * (1 + np.tan(lamA) * np.tan(np.radians(phi)) / Fs))
               ) / np.sum(lamW * np.sin(lamA))
           )

where:

c = 10
phi = 30
lamela1 = {'W':887.36, 'alpha':np.radians(46.71), 'L':19.325}
lamela2 = {'W':1624.8, 'alpha':np.radians(22.054), 'L':14.297}
lamela3 = {'W':737.43, 'alpha':np.radians(1.9096), 'L':13.258}
lamele = [lamela1, lamela2, lamela3]

lamW = np.array([[i['W'] for i in lamele]])
lamA = np.array([[i['alpha'] for i in lamele]])
lamL = np.array([[i['L'] for i in lamele]])

Now, I solved this by making a simple recursion:

def iterf(f, Fsi):
    if np.isclose(Fsi, fn4(Fsi)) == True:
        return Fsi
    else:
        return iterf(f, fn4(Fsi))

which gives

iterf(fn4, 1)
1.8430

But when I try to use scipy.optimize.fsolve (as fsolve(fn4, 1)), then i get:

RuntimeWarning: divide by zero encountered in true_divide
  / (np.cos(lamA) * (1 + np.tan(lamA) * np.tan(np.radians(phi)) / Fs))

I don't know how exactly does fsolve work (which I'd like to learn), but my guess is that for the optimization it needs to input Fs as 0 at one point, which produces the zero-division. How can I use a function like this one to get the output?

EDIT:
The theoretical background is stability of soil slopes. The Fs, which is the result, is the factor of safety - a measure whose values is >1 for stable slopes, and <=1 for unstable slopes. The chosen method is really simple as it requires the solution of only this one equation. The slope is discretized with a number of "lamellae" or "columns", in this case only 3 for simplicity.

Solution

It looks like you're trying to solve fn4(x) == x, that is, find a fixed point of fn4. (If not, please clarify what the recursion is for.) However, the first argument of fsolve is supposed to be the function for which you want a root, so you're actually calling fsolve to solve fn4(x) == 0.

Recommendations:

fsolve is a legacy function; prefer root
root is generalized for multivariable problems, but you have a scalar function, so prefer root_scalar.
Use root_scalar to solve for the fixed point as follows.

import numpy as np

# OP's code
c = 10
phi = 30
lamela1 = {'W':887.36, 'alpha':np.radians(46.71), 'L':19.325}
lamela2 = {'W':1624.8, 'alpha':np.radians(22.054), 'L':14.297}
lamela3 = {'W':737.43, 'alpha':np.radians(1.9096), 'L':13.258}
lamele = [lamela1, lamela2, lamela3]

lamW = np.array([[i['W'] for i in lamele]])
lamA = np.array([[i['alpha'] for i in lamele]])
lamL = np.array([[i['L'] for i in lamele]])

def fn4(Fs):
    return ( np.sum( (lamL * c + lamW * np.tan(np.radians(phi)))
                 / (np.cos(lamA) * (1 + np.tan(lamA) * np.tan(np.radians(phi)) / Fs))
               ) / np.sum(lamW * np.sin(lamA))
           )

# Solution code
from scipy.optimize import fsolve, root, root_scalar
# The `lambda` function returns zero when `fn4(x) == x`
res = root_scalar(lambda x: fn4(x) - x, x0=1)
# or, if you want to restrict the search within a bracket
# res = root_scalar(lambda x: fn4(x) - x, bracket=(0.1, 3))
np.testing.assert_allclose(res.root, fn4(res.root))