Error using symengine>@()0.0 Too many input arguments.

Hello,

I am receiving this error when I run my code:

Error using symengine>@()0.0
Too many input arguments.

Error in odefunction (line 24)
bdip = [s(4); s(5); s(6); (q/m_e)*(Ex(s(1),s(2),s(3)) + s(5)*Bz(s(1),s(2),s(3)) - s(6)*By(s(1),s(2),s(3))); (q/m_e)*(Ey(s(1),s(2),s(3)) +
s(6)*Bx(s(1),s(2),s(3)) - s(4)*Bz(s(1),s(2),s(3))); (q/m_e)*(Ez(s(1),s(2),s(3)) + s(4)*By(s(1),s(2),s(3)) - s(5)*Bx(s(1),s(2),s(3)))];

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode15s (line 150)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Collisions (line 98)
    	[T,S] = ode15s(@odefunction, tspan, icv); %Matlab's ODE solver

This is a more detailed question regarding a previous post I made about a function file. The function file in question previously is what makes me beleive it is the root cause of the issue but I cannot figure out why. I managed to strip my code down from 800 lines down to just about 110 (Not including the plotting section) while still repoducing the same error. I have been attempting to trouble shoot this error for more than a year now and cannot find a solution. I also attached the script files the main code references directly and the plotting is included with the main script file (not shown directly here in this post). I have also already tried using every ode solver Matlab has to offer but the same error keeps popping up.

Reference:

"Collisions.m" = Main script file

"odefunction.m" = ODE function file (refernced by Collisions.m)

"B_test.m" = Bfield function file (refernced by odefunction.m)

"E_test.m" = Efield function file (refernced by odefunction.m)

"Collisions.m" Code:

clear %Clear all variables in the visible workspace
clearvars -global %Clears all global variables not seen in workspace

%Lines 9-23 are for user defined (initial) conditions

%TSTEP MUST BE ABLE TO GO INTO TFIN AN ODD AMOUNT OF TIMES (Because simulation includes tstep @ 0 so it ends up being even)
tstep = 0.0625E-9; %Defining time step
tfin = 20.25E-8; %Defining final time
particlecount = 10; %# of particles to generate energies for energy generation
s1 = 0.0200001; %Specify minimum radius particles may be generated
s2 = 0.0200002; %Specify maximum radius particles may be generated
%Defined space where electrons are generated (in radians)
phi1rad = (pi/180)*90;
phi2rad = (pi/180)*-90;
ien = 0.2; %Initial energy of all electrons in (eV)
T = 290; %Temperature of netural gas

%Natural Constants and Initial Conditions
m_e = 9.11E-31; %Electron mass
epnaut = 8.854187E-12; %Vacuum permitivity constant
k = 1/(4*pi*epnaut); %Constant fraction in terms of VP
initialenergy = (3/2)*k*T;
evtoj = ien*(1.602E-19);

%UI Setup and file saves
savefile = 'test.mat'; %Save file name for radial density function
datasetnum = 'test.mat'; %Name file where initial electron positions, initial electron velocities, collision position, mfp and energy matrix will be saved

tic %Begin timer

%Create zeros matrices to populate later
energymat = zeros(particlecount,1); %Initial energies assigned to each electron
vmagmat = zeros(particlecount,1); %Initial velocity magnitudes assigned to each electron
intspan = [0:tstep:tfin]; %Total time span
[introw,intcol] = size(intspan); %Measures size of matrix to determine how many times the nested for-loop will run
icvelocity = zeros(particlecount,3); %Initial velocities matrix
icposition = zeros(particlecount,3); %Initial positions matrix

%Lines 47-75 are dedicated to generating the initial conditions for the
%all electrons/particles based on the initial conditions provided by the user

%Generates inital energies for each particle
for pe = 1:particlecount
    energymat(pe,1) = evtoj;
end
%Generate random positions for each particle between s1 and s2
for place3 = 1:particlecount
    %In Spherical: Generates random position with boundaries of: s1<r<s2
    rrand = (s2 - s1)*rand() + s1; %Random # generation for component x between s1 and s2
    thetarand = (2*pi)*rand(); %Random # generation for component y between s1 and s2
    phirand = asin((sin(phi2rad) - sin(phi1rad))*rand() + sin(phi1rad)); %Random # generation for component z between s1 and s2
    [xrand,yrand,zrand] = sph2cart(thetarand,phirand,rrand);
    icposition(place3,1:3) = [xrand yrand zrand];
end
%Generates velocity components for each particle
for place2 = 1:particlecount
    vmagin = sqrt(2*energymat(place2,1)/m_e);
    vmagmat(place2,1) = vmagin;
    vx = (icposition(place2,1)/norm(icposition(place2,1:3)))*vmagin;
    vy = (icposition(place2,2)/norm(icposition(place2,1:3)))*vmagin;
    vz = (icposition(place2,3)/norm(icposition(place2,1:3)))*vmagin;
    icvelocity(place2,1:3) = [vx vy vz];
end
%Generate matrix that will be populated by positions and velocities (Each matrix within structure "Wposandvel" should be of size 2*introw-1 by 6)
for place4 = 1:particlecount
    Wposandvel(place4).matrix = zeros(2*introw-1,6);
end
%Generate matrix that will be populated by collisions (Each matrix within structure "collpos" should be of size introw by 3)
for place5 = 1:particlecount
    collpos(place5).matrix = zeros(introw,3);
end

for p = 1:particlecount %particle count

    %Define initial position and velocity vectors for given particle
    spart = icposition(p,1:3);
    vpart = icvelocity(p,1:3);
    
	for t = 0:1:intcol-2 %complete time interval and time step

        [trow,tcol] = size(t);
		%Defining relative position and velocity
		x = spart(1);
		y = spart(2);
		z = spart(3);
		vx = vpart(1);
		vy = vpart(2);
		vz = vpart(3);
		vmag = sqrt(vpart(1).^2 + vpart(2).^2 + vpart(3).^2);
        
		%Coupled differential equation solver for trajectory within current time step
		icv = [x; y; z; vx; vy; vz]; %Initial conditions
        tspan = [intspan(t+1) ((intspan(t+2)-intspan(t+1))/2)+intspan(t+1) intspan(t+2)]; %Time span with given midpoint value to be stored
    	[T,S] = ode15s(@bdipuniodefun, tspan, icv); %Matlab's ODE solver
        Wposandvel(p).matrix((1+2*t):(3+2*t),(1:6)) = S; %Send solution vector to data structure
        
		%Redfine the velocity and position components to reference on next for-loop run
        spart(1) = Wposandvel(p).matrix((3+2*t),1);
		spart(2) = Wposandvel(p).matrix((3+2*t),2);
		spart(3) = Wposandvel(p).matrix((3+2*t),3);
		vpart(1) = Wposandvel(p).matrix((3+2*t),4);
		vpart(2) = Wposandvel(p).matrix((3+2*t),5);
		vpart(3) = Wposandvel(p).matrix((3+2*t),6);
        
    end
end

"odefunction.m" Code:

function bdip = bdipuniodefun(t,s)

%Using SI units
q = -1.60217662E-19;
m_e = 9.11E-31;

persistent Bx By Bz Ex Ey Ez

if isempty(Bx)
    [Bx, By, Bz] = B_test();
end


if isempty(Ex)
    [Ex, Ey, Ez] = E_test();
end

bdip = [s(4); s(5); s(6); (q/m_e)*(Ex(s(1),s(2),s(3)) + s(5)*Bz(s(1),s(2),s(3)) - s(6)*By(s(1),s(2),s(3))); (q/m_e)*(Ey(s(1),s(2),s(3)) + s(6)*Bx(s(1),s(2),s(3)) - s(4)*Bz(s(1),s(2),s(3))); (q/m_e)*(Ez(s(1),s(2),s(3)) + s(4)*By(s(1),s(2),s(3)) - s(5)*Bx(s(1),s(2),s(3)))];

end

If youre wondering why I defined all the equations in a single line (And not like the example shown in Matlab documentation), this is because the ode solver would either not solve the equation set correctly or throw another input argument error.

"B_test.m" Code:

function [Bx, By, Bz] = B_test()

Bfieldstrength = 0.64; %In (Teslas)
magvol = 3.218E-6; %In (m)
mu0 = (4*pi)*10^-7;
magnetization = (Bfieldstrength*magvol)/mu0;

syms x y z
m = [0,0,magnetization];
r = [x, y, z];
B = (mu0/(4*pi))*((dot(m,r)*r*3)/norm(r)^5) - m/norm(r)^3;
Bx = matlabFunction(B(1));
By = matlabFunction(B(2));
Bz = matlabFunction(B(3));

end

"E_test.m" Code:

function [Ex, Ey, Ez] = E_field()

syms x y z
R_s = 0.02;
V = -100;
epnaut = 8.854187E-12;
k = 1/(4*pi*epnaut);
Q = (V*R_s)/k;
r = [x, y, z];
E = (k*Q/norm(r)^3)*r;
Ex = matlabFunction(E(1));
Ey = matlabFunction(E(2));
Ez = matlabFunction(E(3));

end