Когда вы посмотрите, как fsolve определен в модуле scipy, мы увидим:
def fsolve(func, x0, args=(), fprime=None, full_output=0,
col_deriv=0, xtol=1.49012e-8, maxfev=0, band=None,
epsfcn=None, factor=100, diag=None):
"""
Find the roots of a function.
Return the roots of the (non-linear) equations defined by
``func(x) = 0`` given a starting estimate.
Parameters
----------
func : callable ``f(x, *args)``
A function that takes at least one (possibly vector) argument,
and returns a value of the same length.
'''
Таким образом, ваше входное значение для p должно состоять из такого же количества элементов, которое возвращается вашей функцией. Попробуйте, например:
from scipy.optimize import fsolve
import numpy as np
def equations(p):
x1 = p[0]
x2 = p[1]
return x1-8, x2**2 - 64
x = fsolve(equations, np.array([1, 2]))
print(x)
, что дает 8, 8 в качестве ответа.