issue with sympy lambdify

I just want numerically integrate the Lorenz' famous 3 equations I create this code that seems perfectly work:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import sympy as sp
from IPython.display import display
from sympy.interactive import printing
printing.init_printing(use_latex = True)
#%%
def runge_kutta(funzione,Xn,Dt):
    k1 = funzione(Xn) * Dt 
    k2 = funzione(Xn + (1/2)*k1) * Dt
    k3 = funzione(Xn + (1/2)*k2) * Dt
    k4 = funzione(Xn + k3) * Dt
    Xnuovo = Xn + (1/6)*(k1 + 2*k2 + 2*k3 + k4)
    return Xnuovo
#%%
# flux phase of the trajectory
def flusso_traiettoria_singola(condizioni_iniziali, dt, punti_di_tempo, FUNZA):
    dimensionalita = len(condizioni_iniziali)
    X = np.copy(condizioni_iniziali)
    risultati_integrazione = []
    for array_separati in range(dimensionalita):
        risultati_integrazione.append(np.array([condizioni_iniziali[array_separati]]))
    for r in range(1,punti_di_tempo):
        X = runge_kutta(FUNZA,X,dt)
        for array_cui_aggiungere in range(dimensionalita):
            risultati_integrazione[array_cui_aggiungere] = np.append(risultati_integrazione[array_cui_aggiungere], X[array_cui_aggiungere])
    return risultati_integrazione
#%%
sigma = 10
r = 28
b = 8/3
def Lorenz_mio(X):
    x = X[0]
    y = X[1]
    z = X[2]
    
    dxdt = sigma * (y - x)
    dydt = r*x - y - x*z
    dzdt = x*y - b*z
    vettorederivate = np.array([dxdt,dydt,dzdt])

    return vettorederivate
#%%
# INTEGRATION

condizione_iniziale = np.array([1, 1, 1]) # initial condition
dt = 0.01
numero_di_punti_tempo_voluti = 4000
tempo_finale = numero_di_punti_tempo_voluti * dt
t = np.linspace(0,tempo_finale,numero_di_punti_tempo_voluti)
risultati_integrazione_lorenz = flusso_traiettoria_singola(condizione_iniziale, dt, numero_di_punti_tempo_voluti, Lorenz_mio)

# PLOT
fig = plt.figure(figsize=(8, 6))
ax = Axes3D(fig)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.plot_wireframe(risultati_integrazione_lorenz[0], risultati_integrazione_lorenz[1], risultati_integrazione_lorenz[2], color='green', label='ritratto di fase')
plt.savefig()

ok, now I would like do the same thing but using Numpy (eg I can also calculate the fixed points..) and I would not create a function for Lorenz system, as before, but using instead the one I have create with Numpy, my code is:

# NOW SYMPY

x_simbolo, y_simbolo, z_simbolo = sp.symbols('x y z')
sigma_par, r_par, b_par = sp.symbols(r' \sigma r b')

x_dot_deriv_tempo = sigma_par * ( y_simbolo - x_simbolo )
y_dot_deriv_tempo = r_par * x_simbolo - y_simbolo - x_simbolo * z_simbolo
z_dot_deriv_tempo = x_simbolo * y_simbolo - b_par * z_simbolo
sistema_3_eq_lorenz_x_y_z = x_dot_deriv_tempo, y_dot_deriv_tempo, z_dot_deriv_tempo
display(sistema_3_eq_lorenz_x_y_z)

# FIXED POINTS
punti_fissi_sistema_3_eq_lorenz = sp.solve(sistema_3_eq_lorenz_x_y_z, x_simbolo, y_simbolo, z_simbolo, dict=True)
print('i punti fissi delle 3 equazioni sono')
display(punti_fissi_sistema_3_eq_lorenz)
#%%

sympify_integrazione_numerica_sistema_lorenz = sp.sympify( sistema_3_eq_lorenz_x_y_z )
sigma_int_num = 10
r_int_num = 28
b_int_num = 8/3
sympify_integrazione_numerica_sistema_lorenz = sympify_integrazione_numerica_sistema_lorenz.subs([ (sigma_par, sigma_int_num), (r_par, r_int_num), (b_par, b_int_num) ])
lambdifylorenz = sp.lambdify([x_simbolo, y_simbolo, z_simbolo], sympify_integrazione_numerica_sistema_lorenz, "numpy")

risultati_integrazione_lorenz = flusso_traiettoria_singola(condizione_iniziale, dt, numero_di_punti_tempo_voluti, lambdifylorenz)

# PLOT
fig = plt.figure(figsize=(8, 6))
ax = Axes3D(fig)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.plot_wireframe(risultati_integrazione_lorenz[0], risultati_integrazione_lorenz[1], risultati_integrazione_lorenz[2], color='green', label='ritratto di fase')
ax.legend(loc='lower center')
plt.savefig()

Here comes the problem:

lambdifylorenz(condizione_iniziale) gives error

TypeError: <lambda>() missing 2 required positional arguments: '_Dummy_32' and '_Dummy_33'

since I can't use it as I did with

Lorenz_mio(condizione_iniziale)

which correctly gives:

array([ 0. , 26. , -1.66666667])

Of course this will work

lambdifylorenz(condizione_iniziale[0],condizione_iniziale[1],condizione_iniziale[2])

or simply

lambdifylorenz(1,1,1)

but it isn't what I need since it must work with the call

funzione(Xn)

that appears in the functions runge_kutta() and flusso_traiettoria_singola().

I read the Lambdify webpage but I found nothing: there is an example similar to the following

a = np.array([1, 2, 3])
fxy = sp.lambdify([x,y], [x + 1,y**2], 'numpy')

with fxy(a,a) you have [array([2, 3, 4]), array([1, 4, 9])]

Your code was broken in numerous ways; parameter packing was only one of them. Don't set risultati_integrazione to a list; it should be a pre-allocated Numpy array. Don't write an inner for-loop to append. Don't 8/3; this is symbolic math so you should be using a Rational. Don't call Axes3D directly. Don't call plot_wireframe; simply call plot. Don't call sympify; it's easier to work with this as a Matrix since Sympy will coerce the numerical function to an array instead of a tuple. Don't call .legend(); you only have one series.

import typing

import numpy as np
import sympy
import sympy as sp
from matplotlib import pyplot as plt

"""
Sono necessari più commenti in italiano. Il contenuto stampato dovrebbe essere in italiano. 
I nomi delle variabili e delle funzioni dovrebbero essere in inglese.
"""


def runge_kutta(
    function: typing.Callable,  # funzione
    Xn: np.ndarray,
    dt: float,
) -> np.ndarray:
    k1 = dt*function(Xn)
    k2 = dt*function(Xn + k1*0.5)
    k3 = dt*function(Xn + k2*0.5)
    k4 = dt*function(Xn + k3)
    Xnew = Xn + (k1 + 2*(k2 + k3) + k4)/6  # X nuovo
    return Xnew


def trajectory_flux_phase(
    initial_conditions: np.ndarray,  # condizioni iniziali
    dt: float,
    time_points: int,  # punti di tempo
    function: typing.Callable,  # funza
) -> np.ndarray:
    """flusso traiettoria singola"""

    dimensionality = len(initial_conditions)  # dimensionalita
    X = np.copy(initial_conditions)
    integration_results = np.empty((dimensionality, time_points))  # risultati integrazione
    integration_results[:, 0] = initial_conditions

    for r in range(1, time_points):
        X[:] = runge_kutta(function, X, dt)
        integration_results[:, r] = X  # array cui aggiungere

    return integration_results


def plot(integration_result) -> plt.Figure:
    fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_title('Ritratto di Fase')
    ax.plot(*integration_result)
    return fig


def main() -> None:
    x, y, z, sigma, r, b = sp.symbols(r'x y z \sigma r b')  # simbolo

    # derivati temporali
    x_dot = sigma*(y - x)
    y_dot = r*x - y - x*z
    z_dot = x*y - b*z

    # sistema 3 eq lorenz x y z
    lorenz_system = x_dot, y_dot, z_dot
    print(lorenz_system)

    # punti fissi sistema 3 eq Lorenz
    fixed_point_system = sp.solve(lorenz_system, x, y, z, dict=True)
    print('i punti fissi delle 3 equazioni sono')
    print(fixed_point_system)

    # sympify integrazione numerica sistema lorenz
    system_tuple = sp.Matrix(lorenz_system).T.subs([
        (sigma, 10),
        (r, 28),
        (b, sympy.Rational(8, 3)),
    ])
    function = sp.lambdify([x, y, z], system_tuple)

    def vector_function(params: np.ndarray) -> np.ndarray:
        return function(*params).squeeze()

    # risultati integrazione lorenz
    integration_result = trajectory_flux_phase(
        initial_conditions=np.ones(3),  # condizione iniziale
        dt=0.01,
        time_points=4_000,  # numero di punti tempo voluti
        function=vector_function,
    )

    plot(integration_result)
    plt.show()


if __name__ == '__main__':
    main()

(\sigma*(-x + y), r*x - x*z - y, -b*z + x*y)
i punti fissi delle 3 equazioni sono
[{x: 0, y: 0, z: 0}, {x: -sqrt(b*r - b), y: -sqrt(b*(r - 1)), z: r - 1}, {x: sqrt(b*r - b), y: sqrt(b*(r - 1)), z: r - 1}]

More broadly, you should not be re-implementing Runge-Kutta; instead call into e.g. https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html via solve_ivp. This is more accurate, flexible and fast, especially if you don't use Runge-Kutta and use instead e.g. LSODA with an exact Jacobian:

import functools
import time

import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import check_grad


def lorenz_dt(
    t: float, v: np.ndarray, sigma: float, r: float, b: float,
) -> np.ndarray:
    """derivati temporali"""
    x, y, z = v

    return np.array((
        sigma*(y - x),  # xdot
        x*(r - z) - y,  # ydot
        x*y - b*z,      # zdot
    ))


def lorenz_jacobian(
    t: float, v: np.ndarray, sigma: float, r: float, b: float,
) -> np.ndarray:
    x, y, z = v

    return np.array((
        #   /dx    /dy /dz
        (-sigma, sigma,  0),  # df[0]
        (   r-z,    -1, -x),  # df[1]
        (     y,     x, -b),  # df[2]
    ))


def trajectory_flux_phase(
    sigma: float = 10, r: float = 28, b: float = 8/3,
) -> np.ndarray:
    # Queste quantità non sono esatte. L'algoritmo le modificherà in modo intelligente.
    equivalent_dt = 0.01
    equivalent_n = 4_000
    t_end = equivalent_dt*equivalent_n

    initial_value = np.ones(3)  # condizione iniziale
    error = check_grad(
        func=functools.partial(
            lorenz_dt,
            0, sigma=sigma, r=r, b=b,
        ),
        grad=functools.partial(
            lorenz_jacobian,
            0, sigma=sigma, r=r, b=b,
        ),
        x0=initial_value,
    )
    assert error < 1e-6

    t0 = time.perf_counter()
    result = solve_ivp(
        fun=functools.partial(lorenz_dt, sigma=sigma, r=r, b=b),
        jac=functools.partial(lorenz_jacobian, sigma=sigma, r=r, b=b),
        method='LSODA',
        t_span=(0, t_end),
        y0=initial_value,
        rtol=1e-8,
    )
    t1 = time.perf_counter()
    if not result.success:
        raise ValueError(result.message)
    print(result.message)
    print('numero di chiamate di funzione:', result.nfev)
    print("numero di punti temporali scelti dall'algoritmo:", result.y.shape[1])
    print(f'tempo di esecuzione: {t1 - t0:.4f} s')
    print()

    return result.y


def plot(integration_result: np.ndarray) -> plt.Figure:
    fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_title('Ritratto di Fase')
    ax.plot(*integration_result)
    return fig


def main() -> None:
    integration_result = trajectory_flux_phase()
    plot(integration_result)
    plt.show()


if __name__ == '__main__':
    main()