Projectile motion 3d odes

Question

I am trying to solve this set of differential equations with solve_ivp in Python. Winds are horizontal two dimensional vectors, so I need differential equation solved in y dimension.

This is the error I get:

ValueError: not enough values to unpack (expected 6, got 4)

Is there any other way of solving it? Here is my code:(just projectile motion in x-z plane, distance as function of angle). For the example of two such initial angles (with which we hit approximately the same point) I am trying to calculate the dispersion of hits in the wind according to the point hit in no wind. Thanks for any help or pointers

from random import randrange
from statistics import stdev
from turtle import title
import numpy as np
import scipy as sp
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

    
m = 30
v0 = 301
g = 9.81
ro = 1.28
A = 0.0026
c = 1
k = 0.5 * c * ro * A
B = k / m
V = v0
c_max = 10
c_windx = 0
c_windy = 0

def dSdt(t, S, B, c_windx, c_windy):
    x, vx, y, vy, z, vz = S
    speed = np.hypot(np.hypot((vx + c_windx),(vy + c_windy)), vz)
    return [vx, -B * speed * (vx + c_windx), vy, -B * speed * (vy + c_windy), vz, -B * speed * vz - g]

def distance(angle, B, c_windx, V=v0, t=600):
    v0x = V * np.cos(np.radians(angle))
    v0y = 0
    v0z = V * np.sin(np.radians(angle))
    sol = solve_ivp(dSdt, [0, t], y0 =[0,v0x,0,v0y,0,v0z], t_eval=np.linspace(0,t,10000), args=(B, c_windx,c_windy), atol=1e-7, rtol=1e-4)
    justabove = np.where(np.diff(np.sign(sol.y[2])) < 0)[0][0]
    justunder = justabove + 1
    xloc = (sol.y[0][justabove] + sol.y[0][justunder])/2 
    return xloc

angles = np.linspace(0, 90, 200)
xlocs = np.vectorize(distance)(angles, B=B, c_windx = c_windx)


plt.plot(angles, xlocs)
plt.xlabel('angles []')
plt.ylabel('distance [m]') 
plt.axvline(angles[np.argmax(xlocs)], ls='--', color='r')
plt.text(angles[np.argmax(xlocs)] + 0.2 ,0, round(angles[np.argmax(xlocs)], 2), rotation=90, color='r') 
plt.title('Distance(angle)')
plt.show()

alfa = 15 #its optional

def same_distance_alfa(B,c_windx, V=v0, t=600):
        for a in np.arange(alfa + 5, 90, 0.1):            
            if abs(distance(alfa, B, c_windx, V=v0, t=600) - distance(a, B, c_windx, V=v0, t=600) < 10):
                print((distance(15, B, c_windx, V=v0, t=600) - distance(a, B, c_windx, V=v0, t=600)))
                return a


beta = same_distance_alfa(B)