Problem solving simultaneous ODEs using 4th order Runge-Kutta method (Too many input arguments.)

I am trying to solve 6 ODEs simulataneously (three 2nd order and three 1st order). I am solving it using 4th order Runge-Kutta method. But I am getting an error: Too many input arguments.

Please suggest. My code is following.

tic
clear;
close all;
clc;
ae=10 ;     % aspect ratio
ds=1*10^-3 ;  % equivalent dia
ecc=((ae.^2)-1)./((ae.^2)+1) ;   % eccentricity
g=9.81;
nu=0.89*10^-6 ;
A=0.2;   % wave amplitude
k=0.33/A; % k*A <= 0.33 for Linear wave theory to hold correct
w1=sqrt(g*k); % angular frequency
T=2*pi/w1 ;    % time period
B=1.1;
theta_ini= 0  ;  % azimuthal angle matrix 
phi_ini= pi/4 ;    % polar angle matrix 
vp= 0.1*w1*A ;    % Particle intrinsic velocity (m/s)
H=3700;  % depth of ocean
Fr=w1*A./sqrt(g*H);  % Froude number
z0= 3 ;  

% Time span for evolution of dynamics
t0 = 0; %seconds
tf =20*T; %seconds
dt=0.05 ;
t=t0:dt:tf ;
c1=18*nu/(B*ds^2) ;

% mass tensor in body axes
alpha0=(1./((ae.^2)-1).^(3/2)).*(2*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
beta0=(0.5./((ae.^2)-1).^(3/2)).*(2.*(ae.^2).*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
gamma0=beta0 ;
C11d=alpha0./(2-alpha0) ;
C22d=beta0./(2-beta0) ;
C33d=gamma0./(2-gamma0) ;

% resistance tensor in body axes
K11d=((4/3).*(ae.^(-1/3)).*(1-ae.^2))./(ae-((2*(ae.^2)-1).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K22d=((8/3)*(ae.^(-1/3)).*((ae.^2)-1))./(ae+((2*(ae.^2)-3).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K33d=K22d ;

K33_ini=K11d*cos(phi_ini)*cos(phi_ini)+K22d*sin(phi_ini)*sin(phi_ini) ; % K33 @t=0

% X=[x; vx; y; vy; z; vz; px; py; pz];

% RK4
fx = @(x,vx,y,vy,z,vz,px,py,pz) vx; 
fvx = @(x,vx,y,vy,z,vz,px,py,pz) (1./(1+(1/B).*(C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz))))).*((1/B).*((w1^2)*A*exp(k*z)*sin(k*x-w1*t)+...
    w1*A*exp(k*z)*cos(k*x-w1*t).*(-k*w1*A*exp(k*z)*sin(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t)))+...
    (1/B)*((C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz)))*(w1^2)*A*exp(k*z)*sin(k*x-w1*t))-...
    c1*((K11d.*sin(acos(pz)).*sin(acos(pz))+K22d.*cos(acos(pz)).*cos(acos(pz))).*(vx-w1*A*exp(k*z)*cos(k*x-w1*t)))); 
fy = @(x,vx,y,vy,z,vz,px,py,pz) vy; 
fvy = @(x,vx,y,vy,z,vz,px,py,pz) (1./(1+(1/B).*C22d)).*(-c1.*(K22.*vy));
fz = @(x,vx,y,vy,z,vz,px,py,pz) vz;
fvz = @(x,vx,y,vy,z,vz,px,py,pz) (1./(1+(1/B).*(C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz))))).*((1/B).*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)+...
    w1*A*exp(k*z)*cos(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*sin(k*x-w1*t)))+...
    (1/B)*((C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz)))*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)))-...
    c1.*((K11d.*cos(acos(pz)).*cos(acos(pz))+K22d.*sin(acos(pz)).*sin(acos(pz))).*(vz-w1*A*exp(k*z)*sin(k*x-w1*t)))-(1-1/B).*g); 
fpx = @(x,vx,y,vy,z,vz,px,py,pz) ecc.*((-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz)-ecc.*px.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
    pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
fpy = @(x,vx,y,vy,z,vz,px,py,pz) -ecc.*py.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+...
    (px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
fpz = @(x,vx,y,vy,z,vz,px,py,pz) ecc.*((k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz)-ecc.*pz.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
    pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));

h=dt;

x=zeros(length(t),1);
vx=zeros(length(t),1);
px=zeros(length(t),1);
y=zeros(length(t),1);
vy=zeros(length(t),1);
py=zeros(length(t),1);
z=zeros(length(t),1);
vz=zeros(length(t),1);
pz=zeros(length(t),1);

% Initial field velocity @t=0, x=0, z=0
u0x=w1*A;   
u0y=0;
u0z=0; 

% Initial particle velocity 
v_ini_1p=u0x ;
v_ini_2p=u0y ;
v_ini_3p=-((B-1)*g*ds^2)/(18*nu*K33_ini)+u0z ;

% Initial particle orientation
initial_p1 = sin(phi_ini).*cos(theta_ini);
initial_p2 = sin(phi_ini).*sin(theta_ini);
initial_p3 = cos(phi_ini);

x(1) = 0 ;
vx(1) = v_ini_1p ;
y(1) = 0 ;
vy(1) = v_ini_2p ;
z(1) = -z0 ;
vz(1) = v_ini_3p ;
px(1) = initial_p1 ;
py(1) = initial_p2 ;
pz(1) = initial_p3 ;

for i=1: length(t)-1                              
    k_1x = fx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1vx = fvx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1y = fy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1vy = fvy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1z = fz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1vz = fvz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    
    K1y = Fy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    K1z = Fz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1px = fpx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1py = fpy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1pz = fpz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    
    k_2x = fx(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);   
    k_2vx = fvx(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2y = fy(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2vy = fvy(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2z = fz(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2vz = fvz(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2px = fpx(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2py = fpy(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2pz = fpz(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    
    k_3x = fx(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);   
    k_3vx = fvx(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
    k_3y = fy(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);   
    k_3vy = fvy(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
    k_3z = fz(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);   
    k_3vz = fvz(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
    k_3px = fpx(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);   
    k_3py = fpy(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
    k_3pz = fpz(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
    
    k_4x = fx(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
    k_4vx = fvx(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
    k_4y = fy(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);   
    k_4vy = fvy(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
    k_4z = fz(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);   
    k_4vz = fvz(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
    k_4px = fpx(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);   
    k_4py = fpy(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
    k_4pz = fpz(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
    
    x(i+1) = x(i) + (1/6)*(k_1x+2*k_2x+2*k_3x+k_4x)*h;  
    vx(i+1) = vx(i) + (1/6)*(k_1vx+2*k_2vx+2*k_3vx+k_4vx)*h;  
    px(i+1) = px(i) + (1/6)*(k_1px+2*k_2px+2*k_3px+k_4px)*h;  
    
    y(i+1) = y(i) + (1/6)*(k_1y+2*k_2y+2*k_3y+k_4y)*h;  
    vy(i+1) = vy(i) + (1/6)*(k_1vy+2*k_2vy+2*k_3vy+k_4vy)*h;  
    py(i+1) = py(i) + (1/6)*(k_1py+2*k_2py+2*k_3py+k_4py)*h;  
    
    z(i+1) = z(i) + (1/6)*(k_1z+2*k_2z+2*k_3z+k_4z)*h;  
    vz(i+1) = vz(i) + (1/6)*(k_1vz+2*k_2vz+2*k_3vz+k_4vz)*h;  
    pz(i+1) = pz(i) + (1/6)*(k_1pz+2*k_2pz+2*k_3pz+k_4pz)*h;  
    
end

Error using solution>@(x,vx,y,vy,z,vz,px,py,pz)vx
Too many input arguments.

[James Tursa]

I would suggest you rewrite this code using vectors for your states instead of individual variables. E.g., define a state vector y as follows:

y(1) = x

y(2) = y

y(3) = z

y(4) = vx

y(5) = vy

y(6) = vz

y(7) = px

y(8) = py

y(9) = pz

Then define one derivative function handle:

f = @(t,y) = [y(4:6); other stuff ]

Your downstream code will be greatly simplified and easier to debug, and this f can be passed directly into ode45( ) for comparison results.

Answer 1

[Cris LaPierre]

The error is because you define fx with 9 inputs, but call it with 10.

fx = @(x,vx,y,vy,z,vz,px,py,pz) vx;
% vs
k_1x = fx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));

It looks like you left out time when defining fx. Did you mean to do this (add time to function inputs)?

fx = @(t,x,vx,y,vy,z,vz,px,py,pz) vx;

or this (remove time from calling syntax)

k_1x = fx(x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));

It looks like you'll have this same issue with all your functions.

[Ban]

Thanks. Yeah, I missed out the time variable as a typo.

It's too many inputs and too many eqns. I would like to vectorize everything, like this - X=[x,vx,y,vy,z,vx,px,py,pz]; F=@(t,X)[...all the eqns...]; Then vectorize K1,K2,K3,K4.

But the issue is it shows error stating that variables in X are undefined. Please suggest.

[Torsten]

I wonder why you program the numerical method on your own.

MATLAB has built-in efficient and well-tested ODE integrators you can use for your problem.

This way, you can be quite certain that problems you get with the code are due to your equations, not due to the self-programmed Runge-Kutta ODE integrator.

[Cris LaPierre]

Please share the full error message you are seeing (all the red text)

[AkB]

Undefined function or variable 'x'.

Error in del_2 (line 46)

X=[x vx y vy z vz px py pz];

[Walter Roberson]

That line does not exist in your posted code, not even as a comment.

[AkB]

The posted code is detailed. The error "Undefined function or variable 'x'. Error in del_2 (line 46) X=[x vx y vy z vz px py pz];" occurs when I vectorize everything, like this: X=[x,vx,y,vy,z,vx,px,py,pz]; F=@(t, X)[...all the eqns...]; Then vectorize K1,K2,K3,K4.

Answer 2

[Torsten]

I changed

(1./(1+(1/B).*C22d)).*(-c1.*(K22.*vy));

to

(1./(1+(1/B).*C22d)).*(-c1.*(K22d.*vy));

in the definition of fvy.

And you should check these strange cos(acos(...)) expressions in your differential equations.

ae=10 ;     % aspect ratio
ds=1*10^-3 ;  % equivalent dia
ecc=((ae.^2)-1)./((ae.^2)+1) ;   % eccentricity
g=9.81;
nu=0.89*10^-6 ;
A=0.2;   % wave amplitude
k=0.33/A; % k*A <= 0.33 for Linear wave theory to hold correct
w1=sqrt(g*k); % angular frequency
Tp=2*pi/w1 ;    % time period
B=1.1;
theta_ini= 0  ;  % azimuthal angle matrix 
phi_ini= pi/4 ;    % polar angle matrix 
vp= 0.1*w1*A ;    % Particle intrinsic velocity (m/s)
H=3700;  % depth of ocean
Fr=w1*A./sqrt(g*H);  % Froude number
z0= 3 ;  
c1=18*nu/(B*ds^2) ;
% mass tensor in body axes
alpha0=(1./((ae.^2)-1).^(3/2)).*(2*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
beta0=(0.5./((ae.^2)-1).^(3/2)).*(2.*(ae.^2).*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
gamma0=beta0 ;
C11d=alpha0./(2-alpha0) ;
C22d=beta0./(2-beta0) ;
C33d=gamma0./(2-gamma0) ;
% resistance tensor in body axes
K11d=((4/3).*(ae.^(-1/3)).*(1-ae.^2))./(ae-((2*(ae.^2)-1).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K22d=((8/3)*(ae.^(-1/3)).*((ae.^2)-1))./(ae+((2*(ae.^2)-3).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K33d=K22d ;
K33_ini=K11d*cos(phi_ini)*cos(phi_ini)+K22d*sin(phi_ini)*sin(phi_ini) ; % K33 @t=0
% Time span for evolution of dynamics
t0 = 0; %seconds
tf =20*Tp; %seconds
dt=0.05 ;
tspan=t0:dt:tf ;

% Initial field velocity @t=0, x=0, z=0
u0x=w1*A;   
u0y=0;
u0z=0; 
% Initial particle velocity 
v_ini_1p=u0x ;
v_ini_2p=u0y ;
v_ini_3p=-((B-1)*g*ds^2)/(18*nu*K33_ini)+u0z ;
% Initial particle orientation
initial_p1 = sin(phi_ini).*cos(theta_ini);
initial_p2 = sin(phi_ini).*sin(theta_ini);
initial_p3 = cos(phi_ini);
x = 0 ;
vx = v_ini_1p ;
y = 0 ;
vy = v_ini_2p ;
z = -z0 ;
vz = v_ini_3p ;
px = initial_p1 ;
py = initial_p2 ;
pz = initial_p3 ;
Y0 = [x;vx;y;vy;z;vz;px;py;pz];
[T,Y] = ode45(@fun,tspan,Y0);
plot(T,Y(:,1))

function dY = fun(t,Y)
% Y=[x; vx; y; vy; z; vz; px; py; pz];
x = Y(1);
vx = Y(2);
y = Y(3);
vy = Y(4);
z = Y(5);
vz = Y(6);
px = Y(7);
py = Y(8);
pz = Y(9);

ae=10 ;     % aspect ratio
ds=1*10^-3 ;  % equivalent dia
ecc=((ae.^2)-1)./((ae.^2)+1) ;   % eccentricity
g=9.81;
nu=0.89*10^-6 ;
A=0.2;   % wave amplitude
k=0.33/A; % k*A <= 0.33 for Linear wave theory to hold correct
w1=sqrt(g*k); % angular frequency
T=2*pi/w1 ;    % time period
B=1.1;
theta_ini= 0  ;  % azimuthal angle matrix 
phi_ini= pi/4 ;    % polar angle matrix 
vp= 0.1*w1*A ;    % Particle intrinsic velocity (m/s)
H=3700;  % depth of ocean
Fr=w1*A./sqrt(g*H);  % Froude number
z0= 3 ;  
c1=18*nu/(B*ds^2) ;
% mass tensor in body axes
alpha0=(1./((ae.^2)-1).^(3/2)).*(2*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
beta0=(0.5./((ae.^2)-1).^(3/2)).*(2.*(ae.^2).*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
gamma0=beta0 ;
C11d=alpha0./(2-alpha0) ;
C22d=beta0./(2-beta0) ;
C33d=gamma0./(2-gamma0) ;
% resistance tensor in body axes
K11d=((4/3).*(ae.^(-1/3)).*(1-ae.^2))./(ae-((2*(ae.^2)-1).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K22d=((8/3)*(ae.^(-1/3)).*((ae.^2)-1))./(ae+((2*(ae.^2)-3).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K33d=K22d ;
K33_ini=K11d*cos(phi_ini)*cos(phi_ini)+K22d*sin(phi_ini)*sin(phi_ini) ; % K33 @t=0

dY(1) =  vx; 
dY(2) =  (1./(1+(1/B).*(C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz))))).*((1/B).*((w1^2)*A*exp(k*z)*sin(k*x-w1*t)+...
         w1*A*exp(k*z)*cos(k*x-w1*t).*(-k*w1*A*exp(k*z)*sin(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t)))+...
         (1/B)*((C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz)))*(w1^2)*A*exp(k*z)*sin(k*x-w1*t))-...
         c1*((K11d.*sin(acos(pz)).*sin(acos(pz))+K22d.*cos(acos(pz)).*cos(acos(pz))).*(vx-w1*A*exp(k*z)*cos(k*x-w1*t)))); 
dY(3) = vy; 
dY(4) =  (1./(1+(1/B).*C22d)).*(-c1.*(K22d.*vy));
dY(5) = vz;
dY(6) =  (1./(1+(1/B).*(C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz))))).*((1/B).*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)+...
         w1*A*exp(k*z)*cos(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*sin(k*x-w1*t)))+...
         (1/B)*((C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz)))*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)))-...
         c1.*((K11d.*cos(acos(pz)).*cos(acos(pz))+K22d.*sin(acos(pz)).*sin(acos(pz))).*(vz-w1*A*exp(k*z)*sin(k*x-w1*t)))-(1-1/B).*g); 
dY(7) =  ecc.*((-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz)-ecc.*px.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
         pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
dY(8) = -ecc.*py.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+...
       (px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
dY(9) = ecc.*((k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz)-ecc.*pz.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
        pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
dY = dY.';
end

[Torsten]

Your code works, too:

ae=10 ;     % aspect ratio
ds=1*10^-3 ;  % equivalent dia
ecc=((ae.^2)-1)./((ae.^2)+1) ;   % eccentricity
g=9.81;
nu=0.89*10^-6 ;
A=0.2;   % wave amplitude
k=0.33/A; % k*A <= 0.33 for Linear wave theory to hold correct
w1=sqrt(g*k); % angular frequency
T=2*pi/w1 ;    % time period
B=1.1;
theta_ini= 0  ;  % azimuthal angle matrix 
phi_ini= pi/4 ;    % polar angle matrix 
vp= 0.1*w1*A ;    % Particle intrinsic velocity (m/s)
H=3700;  % depth of ocean
Fr=w1*A./sqrt(g*H);  % Froude number
z0= 3 ;  

% Time span for evolution of dynamics
t0 = 0; %seconds
tf =20*T; %seconds
dt=0.05 ;
t=t0:dt:tf ;
c1=18*nu/(B*ds^2) ;

% mass tensor in body axes
alpha0=(1./((ae.^2)-1).^(3/2)).*(2*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
beta0=(0.5./((ae.^2)-1).^(3/2)).*(2.*(ae.^2).*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
gamma0=beta0 ;
C11d=alpha0./(2-alpha0) ;
C22d=beta0./(2-beta0) ;
C33d=gamma0./(2-gamma0) ;

% resistance tensor in body axes
K11d=((4/3).*(ae.^(-1/3)).*(1-ae.^2))./(ae-((2*(ae.^2)-1).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K22d=((8/3)*(ae.^(-1/3)).*((ae.^2)-1))./(ae+((2*(ae.^2)-3).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K33d=K22d ;

K33_ini=K11d*cos(phi_ini)*cos(phi_ini)+K22d*sin(phi_ini)*sin(phi_ini) ; % K33 @t=0

% X=[x; vx; y; vy; z; vz; px; py; pz];

% RK4
fx = @(t,x,vx,y,vy,z,vz,px,py,pz) vx; 
fvx = @(t,x,vx,y,vy,z,vz,px,py,pz) (1./(1+(1/B).*(C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz))))).*((1/B).*((w1^2)*A*exp(k*z)*sin(k*x-w1*t)+...
    w1*A*exp(k*z)*cos(k*x-w1*t).*(-k*w1*A*exp(k*z)*sin(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t)))+...
    (1/B)*((C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz)))*(w1^2)*A*exp(k*z)*sin(k*x-w1*t))-...
    c1*((K11d.*sin(acos(pz)).*sin(acos(pz))+K22d.*cos(acos(pz)).*cos(acos(pz))).*(vx-w1*A*exp(k*z)*cos(k*x-w1*t)))); 
fy = @(t,x,vx,y,vy,z,vz,px,py,pz) vy; 
fvy = @(t,x,vx,y,vy,z,vz,px,py,pz) (1./(1+(1/B).*C22d)).*(-c1.*(K22d.*vy));
fz = @(t,x,vx,y,vy,z,vz,px,py,pz) vz;
fvz = @(t,x,vx,y,vy,z,vz,px,py,pz) (1./(1+(1/B).*(C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz))))).*((1/B).*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)+...
    w1*A*exp(k*z)*cos(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*sin(k*x-w1*t)))+...
    (1/B)*((C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz)))*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)))-...
    c1.*((K11d.*cos(acos(pz)).*cos(acos(pz))+K22d.*sin(acos(pz)).*sin(acos(pz))).*(vz-w1*A*exp(k*z)*sin(k*x-w1*t)))-(1-1/B).*g); 
fpx = @(t,x,vx,y,vy,z,vz,px,py,pz) ecc.*((-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz)-ecc.*px.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
    pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
fpy = @(t,x,vx,y,vy,z,vz,px,py,pz) -ecc.*py.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+...
    (px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
fpz = @(t,x,vx,y,vy,z,vz,px,py,pz) ecc.*((k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz)-ecc.*pz.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
    pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));

h=dt;

x=zeros(length(t),1);
vx=zeros(length(t),1);
px=zeros(length(t),1);
y=zeros(length(t),1);
vy=zeros(length(t),1);
py=zeros(length(t),1);
z=zeros(length(t),1);
vz=zeros(length(t),1);
pz=zeros(length(t),1);

% Initial field velocity @t=0, x=0, z=0
u0x=w1*A;   
u0y=0;
u0z=0; 

% Initial particle velocity 
v_ini_1p=u0x ;
v_ini_2p=u0y ;
v_ini_3p=-((B-1)*g*ds^2)/(18*nu*K33_ini)+u0z ;

% Initial particle orientation
initial_p1 = sin(phi_ini).*cos(theta_ini);
initial_p2 = sin(phi_ini).*sin(theta_ini);
initial_p3 = cos(phi_ini);

x(1) = 0 ;
vx(1) = v_ini_1p ;
y(1) = 0 ;
vy(1) = v_ini_2p ;
z(1) = -z0 ;
vz(1) = v_ini_3p ;
px(1) = initial_p1 ;
py(1) = initial_p2 ;
pz(1) = initial_p3 ;

for i=1: length(t)-1                              
    k_1x = fx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1vx = fvx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1y = fy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1vy = fvy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1z = fz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1vz = fvz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    
    K1y = fy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    K1z = fz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1px = fpx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1py = fpy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    k_1pz = fpz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
    
    k_2x = fx(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);   
    k_2vx = fvx(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2y = fy(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2vy = fvy(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2z = fz(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2vz = fvz(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2px = fpx(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2py = fpy(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    k_2pz = fpz(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
    
    k_3x = fx(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);   
    k_3vx = fvx(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
    k_3y = fy(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);   
    k_3vy = fvy(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
    k_3z = fz(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);   
    k_3vz = fvz(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
    k_3px = fpx(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);   
    k_3py = fpy(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
    k_3pz = fpz(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
    
    k_4x = fx(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
    k_4vx = fvx(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
    k_4y = fy(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);   
    k_4vy = fvy(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
    k_4z = fz(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);   
    k_4vz = fvz(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
    k_4px = fpx(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);   
    k_4py = fpy(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
    k_4pz = fpz(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
    
    x(i+1) = x(i) + (1/6)*(k_1x+2*k_2x+2*k_3x+k_4x)*h;  
    vx(i+1) = vx(i) + (1/6)*(k_1vx+2*k_2vx+2*k_3vx+k_4vx)*h;  
    px(i+1) = px(i) + (1/6)*(k_1px+2*k_2px+2*k_3px+k_4px)*h;  
    
    y(i+1) = y(i) + (1/6)*(k_1y+2*k_2y+2*k_3y+k_4y)*h;  
    vy(i+1) = vy(i) + (1/6)*(k_1vy+2*k_2vy+2*k_3vy+k_4vy)*h;  
    py(i+1) = py(i) + (1/6)*(k_1py+2*k_2py+2*k_3py+k_4py)*h;  
    
    z(i+1) = z(i) + (1/6)*(k_1z+2*k_2z+2*k_3z+k_4z)*h;  
    vz(i+1) = vz(i) + (1/6)*(k_1vz+2*k_2vz+2*k_3vz+k_4vz)*h;  
    pz(i+1) = pz(i) + (1/6)*(k_1pz+2*k_2pz+2*k_3pz+k_4pz)*h;  
    
end
plot(t,x)

[AkB]

Just a query: when I used the same eqns with another different code, it gives different results. Please check.

tic
clear;
close all;
clc;
ae=10 ;     % aspect ratio
ds=1*10^-3 ;  % equivalent dia
ecc=((ae.^2)-1)./((ae.^2)+1) ;   % eccentricity
g=9.81;
nu=0.89*10^-6 ;
A=0.2;   % wave amplitude
k=0.33/A; % k*A <= 0.33 for Linear wave theory to hold correct
w1=sqrt(g*k); % angular frequency
T=2*pi/w1 ;    % time period
B=1.1;
theta_ini= 0  ;  % azimuthal angle matrix 
phi_ini= pi/4 ;    % polar angle matrix 
vp= 0.1*w1*A ;    % Particle intrinsic velocity (m/s)


theta_ini_matrix =[0 pi/4 pi/2] ;  % azimuthal angle matrix 
phi_ini_matrix =[0 pi/4 pi/2] ;    % polar angle matrix 
vp_matrix=[10^-1 10^-2 10^-3].*w1*A ;    % Particle intrinsic velocity (m/s)

H=3700;  % depth of ocean
Fr=w1*A./sqrt(g*H);  % Froude number

% Initial depth of particles
z0= 3 ;  

% time scale of photoaxis
tau_s_matrix =[0.1 1 10] ; 

% Time span for evolution of dynamics
t0 = 0; %seconds
tf =20*T; %seconds
interval = [t0 tf];  %time interval  
n_points=200;
tValues = linspace(interval(1),interval(2),tf*n_points);

syms p1(t) p2(t) p3(t) x(t) y(t) z(t)  

% mass tensor in body axes
alpha0=(1./((ae.^2)-1).^(3/2)).*(2*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
beta0=(0.5./((ae.^2)-1).^(3/2)).*(2.*(ae.^2).*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
gamma0=beta0 ;
C11d=alpha0./(2-alpha0) ;
C22d=beta0./(2-beta0) ;
C33d=gamma0./(2-gamma0) ;

% resistance tensor in body axes
K11d=((4/3).*(ae.^(-1/3)).*(1-ae.^2))./(ae-((2*(ae.^2)-1).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K22d=((8/3)*(ae.^(-1/3)).*((ae.^2)-1))./(ae+((2*(ae.^2)-3).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K33d=K22d ;

K33_ini=K11d*cos(phi_ini)*cos(phi_ini)+K22d*sin(phi_ini)*sin(phi_ini) ; % K33 @t=0

% Defining the governing equations 
c1=18*nu/(B*ds^2) ;

vx=diff(x,t);
vy=diff(y,t);
vz=diff(z,t);

% Maxey–Riley Equation (passive)
eq1 = diff(x,t,t) == (1./(1+(1/B).*(C11d.*sin(acos(p3)).*sin(acos(p3))+C22d.*cos(acos(p3)).*cos(acos(p3))))).*((1/B).*((w1^2)*A*exp(k*z)*sin(k*x-w1*t)+...
    w1*A*exp(k*z)*cos(k*x-w1*t).*(-k*w1*A*exp(k*z)*sin(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t)))+...
    (1/B)*((C11d.*sin(acos(p3)).*sin(acos(p3))+C22d.*cos(acos(p3)).*cos(acos(p3)))*(w1^2)*A*exp(k*z)*sin(k*x-w1*t))-...
    c1*((K11d.*sin(acos(p3)).*sin(acos(p3))+K22d.*cos(acos(p3)).*cos(acos(p3))).*(vx-w1*A*exp(k*z)*cos(k*x-w1*t))));
eq2 = diff(y,t,t) == (1./(1+(1/B).*C22d)).*(-c1.*(K22d.*vy));
eq3 = diff(z,t,t) == (1./(1+(1/B).*(C11d.*cos(acos(p3)).*cos(acos(p3))+C22d.*sin(acos(p3)).*sin(acos(p3))))).*((1/B).*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)+...
    w1*A*exp(k*z)*cos(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*sin(k*x-w1*t)))+...
    (1/B)*((C11d.*cos(acos(p3)).*cos(acos(p3))+C22d.*sin(acos(p3)).*sin(acos(p3)))*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)))-...
    c1.*((K11d.*cos(acos(p3)).*cos(acos(p3))+K22d.*sin(acos(p3)).*sin(acos(p3))).*(vz-w1*A*exp(k*z)*sin(k*x-w1*t)))-(1-1/B).*g); 

% Jeffery’s Equation 
eq4 = diff(p1,t) == ecc.*((-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p1+(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p3)-ecc.*p1.*((p1.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p1+...
    p3.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p1)+(p1.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p3+p3.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p3));
eq5 = diff(p2,t) == -ecc.*p2.*((p1.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p1+p3.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p1)+...
    (p1.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p3+p3.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p3));
eq6 = diff(p3,t) == ecc.*((k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p1+(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p3)-ecc.*p3.*((p1.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p1+...
    p3.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p1)+(p1.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p3+p3.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p3));


vars = [x(t); y(t); z(t); p1(t); p2(t); p3(t)];
V_p = odeToVectorField([eq2 eq1 eq3 eq4 eq5 eq6])

V_p = 

M_p = matlabFunction(V_p,'vars', {'t','Y'});


% Initial field velocity @t=0, x=0, z=0
u0x=w1*A;   
u0y=0;
u0z=0; 

% Initial particle velocity 
v_ini_1p=u0x ;
v_ini_2p=u0y ;
v_ini_3p=-((B-1)*g*ds^2)/(18*nu*K33_ini)+u0z ;

% Initial particle orientation
initial_p1 = sin(phi_ini).*cos(theta_ini);
initial_p2 = sin(phi_ini).*sin(theta_ini);
initial_p3 = cos(phi_ini);

%# Solving for passive particle 
y0_p=[0 v_ini_1p 0 v_ini_2p 0 v_ini_3p initial_p1 initial_p2 initial_p3]; %Initial Conditions for passive particle 
ySol_p = ode45(M_p,interval,y0_p); %solution for passive particle 

x_passive = deval(ySol_p,tValues,1);
y_passive = deval(ySol_p,tValues,3);
z_passive = -z0+deval(ySol_p,tValues,5);

v1_passive = deval(ySol_p,tValues,2);
v2_passive = deval(ySol_p,tValues,4);
v3_passive = deval(ySol_p,tValues,6);

p1_passive = deval(ySol_p,tValues,7);
p2_passive = deval(ySol_p,tValues,8);
p3_passive = deval(ySol_p,tValues,9);

phi_passive=acos(p3_passive) ;
theta_passive=asin(real(p2_passive./sin(phi_passive)));


%# plotting of results
plot(tValues,x_passive,'LineWidth',2)  
xlabel('$t$','FontSize',20,'FontWeight','bold', 'Interpreter','latex');
ylabel('$x$', 'FontSize',20,'FontWeight','bold', 'Interpreter','latex');
set(gca,'FontSize',15);

toc

Elapsed time is 15.046097 seconds.

[Torsten]

Well, I won't go to search now where the equations differ. But they must differ - that's for sure.

Both systems are solved using ODE45. If the equations were identical, ODE45 would have produced the same solution.