def METHOD_NAME(AA, timeseries, N_max, sign, symNx = 2, throw_out_modes=False):
    """ Constructs the matrix A and the vector b from a timeseries of toy
    action-angles AA to solve for the vector x = (theta_0, theta_1, theta_2, omega_1,
    omega_2, omega_3, dSdx..., dSdy..., dSdz...) where x contains all derivatives
    of the Fourier components of the generating function with |n| < N_max """
    # First unroll angles
    angs = unroll_angles(AA.T[3:].T, sign)
    # Same considerations as above
    symNz = 2
    NNx = range(-N_max, N_max+1, symNx)
    NNy = range(-N_max, N_max+1, symNz)
    NNz = range(-N_max, N_max+1, symNz)
    n_vectors = np.array([[i, j, k] for (i, j, k) in product(NNx, NNy, NNz)
                          if(not(i==0 and j==0 and k==0)    # exclude zero vector
                             and (k>0                          # northern hemisphere
                                  or (k==0 and j>0)                 # half of x-y plane
                                  or (k==0 and j==0 and i>0))       # half of x axis
                             and np.sqrt(i*i+j*j+k*k)<=N_max     # inside sphere
                             )])
    if(throw_out_modes):
        n_vectors = np.delete(n_vectors, check_each_direction(n_vectors, angs), axis=0)
    nv = len(n_vectors)
    n = 3*nv+6
    b = np.zeros(shape=(n, ))
    a = np.zeros(shape=(n, n))
    a[:3, :3]=len(AA)*np.identity(3)
    a[:3, 3:6]=np.sum(timeseries)*np.identity(3)
    a[3:6, :3]=a[:3, 3:6]
    a[3:6, 3:6]=np.sum(timeseries*timeseries)*np.identity(3)
    for i, j in zip(angs, timeseries):
        a[6:6+nv, 0]+=-2.*np.sin(np.dot(n_vectors, i))
        a[6:6+nv, 3]+=-2.*j*np.sin(np.dot(n_vectors, i))
        a[6:6+nv, 6:6+nv]+=4.*np.outer(np.sin(np.dot(n_vectors, i)), np.sin(np.dot(n_vectors, i)))
        b[:3]+=i
        b[3:6]+=j*i
        b[6:6+nv]+=-2.*i[0]*np.sin(np.dot(n_vectors, i))
        b[6+nv:6+2*nv]+=-2.*i[1]*np.sin(np.dot(n_vectors, i))
        b[6+2*nv:6+3*nv]+=-2.*i[2]*np.sin(np.dot(n_vectors, i))
    a[6+nv:6+2*nv, 1]=a[6:6+nv, 0]
    a[6+2*nv:6+3*nv, 2]=a[6:6+nv, 0]
    a[6+nv:6+2*nv, 4]=a[6:6+nv, 3]
    a[6+2*nv:6+3*nv, 5]=a[6:6+nv, 3]
    a[6+nv:6+2*nv, 6+nv:6+2*nv]=a[6:6+nv, 6:6+nv]
    a[6+2*nv:6+3*nv, 6+2*nv:6+3*nv]=a[6:6+nv, 6:6+nv]
    a[:6, :]=a[:, :6].T
    return np.array(solve(a, b))