Plotting implicit function using fsolve

I am trying to plot an implicit function using scipy.fsolve but cant seem to get it to work.
I want to plot the function z(x,y) where x + y + z + sin(z) = 0 using scipy.fsolve.

import numpy as np
import scipy
import matplotlib.pyplot as plt
f = lambda x,y: scipy.optimize.fsolve(lambda z: x + y + z + np.sin(z), x0 = 0).tolist()[0]
x = np.linspace(-1,1,50)
y = np.linspace(-1,1,40)
X,Y = np.meshgrid(x,y) 
Z = f(X,Y)
fig =plt.figure()
plt.show()

`
This gives me "_minpack.error: Result from function call is not a proper array of floats".

I dont see why this would be since I can call the function on individual values and its return type is float. It seems like the controversial line is Z = f(X,Y).

Thankful for any help!

EDIT
The original error is fixed for the most part however

def z_func(x, y): z_solve = scipy.optimize.fsolve(lambda z: x + y + z + np.sin(z), x0=np.zeros_like(x))
return z_solve.tolist()[0] 
x = np.linspace(-1,1,50) 
y = np.linspace(-1,1,40) 
X,Y = np.meshgrid(x,y) 
X = X.reshape(-1)
Y = Y.reshape(-1) 
Z = z_func(X,Y)
fig =plt.figure()
ax = plt.axes(projection='3d')  
ax.plot_wireframe(X, Y, Z, color = 'red')
plt.show()

now gives me AttributeError: 'float' object has no attribute 'ndim'

> ` This gives me "_minpack.error: Result from function call is not a proper array of floats".

This is because you're sending 2D arrays (X and Y) in the fsolve function. Try reshaping X and Y before the f() call.

X = X.reshape(-1)
Y = Y.reshape(-1)

You'll come across a second error because the x0 is set to a unique 0 value while it should have the same shape as X and Y. Here is an example (you can flat it as a lambda function if you want):

def z_func(x, y):
    z_solve = scipy.optimize.fsolve(lambda z: x + y + z + np.sin(z), x0=np.zeros_like(x))
    return z_solve.tolist()[0]

EDIT following the additional information provided in the comments section

To plot a 3D wireframe, you will X, Y and Z inputs to be 2D arrays. So the steps would be:

• keep the 2D shape of either X or Y (same ones) outputted by the meshgrid function,
• flatten X and Y to get all positions and compute Z for each position.
• reshape X, Y and Z to the initial 2D shape
• call ax.plot_wireframe(X, Y, Z)

# [...]
X, Y = np.meshgrid(x,y) 
grid_shape = X.shape  # same has Y.shape
# Flatten X and Y to get all possible positions
X = X.reshape(-1) 
Y = Y.reshape(-1) 
positions = np.dstack([X, Y])[0]
# Compute Z values for each position
Z = np.array([f(xpos, ypos) for xpos, ypos in positions])
# Reshape everything to be plotted in 3D axis
X = X.reshape(grid_shape)
Y = Y.reshape(grid_shape)
Z = Z.reshape(grid_shape)
fig = plt.figure() 
ax = plt.axes(projection='3d') 
ax.plot_wireframe(X, Y, Z, color = 'red')

You should get the following figure :