Matlab simulating winnowing using Euler and RK4

Question

I am trying to make a winnowing simulator using matlab. Here is my code.

My Euler simulation works as expected but I am unable to figure out the error in rk4 code. Here is my full code for reference. My rk4 output doesnt show grain trajectory and chaff is in freefall. Absolute noobie with matlab so please help. I am unfamiliar with most of the concepts and trying to learn so sometimes I miss out understanding. So even pointers to where I am missing stuff would be great.

% Simulation parameters
t_start = 0;
dt = 0.01;
t_end = 10;

time = t_start:dt:t_end;

v_0 = 0;
w_0 = 0;
global u_0;
u_0 = 0.5; % Increased horizontal force
x_0 = 0.5;
y_0 = 0.5;
x_c = 0.55;
y_c = -0.5;

% Allocate arrays
y_pos = zeros(length(time));
v_vel = zeros(length(time));
x_pos = zeros(length(time));
w_vel = zeros(length(time));

% Initial conditions
y_pos(1) = y_0; 
v_vel(1) = v_0; 
x_pos(1) = x_0; 
w_vel(1) = 0;

% Constants and parameters
h_slot = 0.1; 
density_air = 1.2; 
density_particle_grain = 750;
viscosity_air = 1.8 * 10^-5;
diameter_particle_grain = 0.0025; 
density_particle_chaff = 50; 
diameter_particle_chaff = 0.00325;

gravity = -9.81; velocity_final = 0;

volume_particle_grain = ((pi/6) * diameter_particle_grain^3);
volume_particle_chaff = ((pi/6) * diameter_particle_chaff^3);

% Compute grain trajectory using Euler's method
for i = 1:numel(time)-1
    if y_pos(i) < y_c
        break; % Stop loop when y_pos(i) < y_c
    end
    
    % Equations for grain
    force_gravity_grain = density_particle_grain * volume_particle_grain * gravity;
    force_buoyancy_grain = (-density_air * volume_particle_grain * gravity);
    mass_grain = density_particle_grain * volume_particle_grain;
    fluid_velocity = 6.2 * u_0 * exp(-50 * ((y_pos(i) / x_pos(i))^2)) * sqrt(h_slot / x_pos(i));

    Reynolds_number_grain = density_air * sqrt((fluid_velocity - w_vel(i))^2 + (velocity_final - v_vel(i))^2) * diameter_particle_grain / viscosity_air;
    if Reynolds_number_grain < 800
        drag_coefficient = (24 / Reynolds_number_grain) * (1 + (0.15 * (Reynolds_number_grain)^(0.687)));
    else 
        drag_coefficient = 0.44;
    end 
    drag_force_x = (0.5 * pi * density_air * (diameter_particle_grain^2) * drag_coefficient * sqrt((fluid_velocity - w_vel(i))^2 + (velocity_final - v_vel(i))^2) * (fluid_velocity - w_vel(i)));
    drag_force_y = (0.5 * pi * density_air * (diameter_particle_grain^2) * drag_coefficient * sqrt((fluid_velocity - w_vel(i))^2 + (velocity_final - v_vel(i))^2) * (velocity_final - v_vel(i)));

    force_tangential = drag_force_x / mass_grain;
    force_vertical = (force_gravity_grain + force_buoyancy_grain + drag_force_y) / mass_grain;

    y_pos(i+1) = y_pos(i) + dt * v_vel(i);
    x_pos(i+1) = x_pos(i) + dt * w_vel(i);
    v_vel(i+1) = v_vel(i) + dt * force_vertical;
    w_vel(i+1) = w_vel(i) + dt * force_tangential;
end

xg_euler = x_pos(1:i); % Grain x positions for Euler
yg_euler = y_pos(1:i); % Grain y positions for Euler

% Compute chaff trajectory using Euler's method
y_pos = zeros(length(time));
v_vel = zeros(length(time));
x_pos = zeros(length(time));
w_vel = zeros(length(time));

y_pos(1) = y_0;
v_vel(1) = v_0; 
x_pos(1) = x_0; 
w_vel(1) = 0;

for i = 1:numel(time)-1
    if y_pos(i) < y_c
        break; % Stop loop when y_pos(i) < y_c
    end
    
    % Equations for chaff
    force_gravity_chaff = density_particle_chaff * volume_particle_chaff * gravity;
    force_buoyancy_chaff = (-density_air * volume_particle_chaff * gravity);
    mass_chaff = density_particle_chaff * volume_particle_chaff;
    fluid_velocity = 6.2 * u_0 * exp(-50 * ((y_pos(i) / x_pos(i))^2)) * sqrt(h_slot / x_pos(i));

    Reynolds_number_chaff = density_air * sqrt((fluid_velocity - w_vel(i))^2 + (velocity_final - v_vel(i))^2) * diameter_particle_chaff / viscosity_air;
    if Reynolds_number_chaff < 800
        drag_coefficient = (24 / Reynolds_number_chaff) * (1 + (0.15 * (Reynolds_number_chaff)^(0.687)));
    else 
        drag_coefficient = 0.44;
    end 
    drag_force_x = (0.5 * pi * density_air * (diameter_particle_chaff^2) * drag_coefficient * sqrt((fluid_velocity - w_vel(i))^2 + (velocity_final - v_vel(i))^2) * (fluid_velocity - w_vel(i)));
    drag_force_y = (0.5 * pi * density_air * (diameter_particle_chaff^2) * drag_coefficient * sqrt((fluid_velocity - w_vel(i))^2 + (velocity_final - v_vel(i))^2) * (velocity_final - v_vel(i)));

    force_tangential = drag_force_x / mass_chaff;
    force_vertical = (force_gravity_chaff + force_buoyancy_chaff + drag_force_y) / mass_chaff;

    y_pos(i+1) = y_pos(i) + dt * v_vel(i);
    x_pos(i+1) = x_pos(i) + dt * w_vel(i);
    v_vel(i+1) = v_vel(i) + dt * force_vertical;
    w_vel(i+1) = w_vel(i) + dt * force_tangential;
end

xch_euler = x_pos(1:i); % Chaff x positions for Euler
ych_euler = y_pos(1:i); % Chaff y positions for Euler

% Compute grain trajectory using RK4 method
[yg_rk4, xg_rk4] = rungeKutta4(@computeTrajectory, time, [y_0, v_0, x_0, w_0], y_c, density_air, u_0, h_slot, density_particle_grain, diameter_particle_grain, viscosity_air, gravity, volume_particle_grain);
% Compute chaff trajectory using RK4 method
[ych_rk4, xch_rk4] = rungeKutta4(@computeTrajectory, time, [y_0, v_0, x_0, w_0], y_c, density_air, u_0, h_slot, density_particle_chaff, diameter_particle_chaff, viscosity_air, gravity, volume_particle_chaff);

% Plotting
figure;

% Euler Method
subplot(1, 2, 1);
plot(xg_euler, yg_euler, 'r-', 'LineWidth', 1.5);
hold on;
plot(xch_euler, ych_euler, 'g--', 'LineWidth', 1.5);
scatter(x_c + 0.05, y_c, 50, 'b', 'filled'); % Separation point
plot([x_c - 0.05, x_c + 0.05], [y_c, y_c], 'k--', 'LineWidth', 1.5); % Denote Bin1 and Bin2 region with black dashed line
text(x_c - 0.075, y_c - 0.05, 'Bin1', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'top');
text(x_c + 0.075, y_c - 0.05, 'Bin2', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'top');
xlabel('x (m)');
ylabel('y (m)');
title('Euler Method');
legend('Grain', 'Chaff', 'Separation Point', 'Location', 'Best');
axis equal;
grid on;
hold off;

% RK4 Method
subplot(1, 2, 2);
plot(xg_rk4, yg_rk4, 'r-', 'LineWidth', 1.5);
hold on;
plot(xch_rk4, ych_rk4, 'g--', 'LineWidth', 1.5);
scatter(x_c + 0.05, y_c, 50, 'b', 'filled'); % Separation point
plot([x_c - 0.05, x_c + 0.05], [y_c, y_c], 'k--', 'LineWidth', 1.5); % Denote Bin1 and Bin2 region with black dashed line
text(x_c - 0.075, y_c - 0.05, 'Bin1', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'top');
text(x_c + 0.075, y_c - 0.05, 'Bin2', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'top');
xlabel('x (m)');
ylabel('y (m)');
title('RK4 Method');
legend('Grain', 'Chaff', 'Separation Point', 'Location', 'Best');
axis equal;
grid on;
hold off;

% Function defining the particle motion equations for RK4
function dydt = computeTrajectory(t, y, yc, denf, u0, h, denp, dp, visf, g, Vp)
    % Extracting variables
    y_particle = y(1);
    v_particle = y(2);
    x_particle = y(3);
    w_particle = y(4);

    % Equations of motion
    Fg = denp * Vp * g;
    Fb = -denf * Vp * g;
    m_particle = denp * Vp;
    uf = 6.2 * u0 * exp(-50 * ((y_particle / x_particle)^2)) * sqrt(h / x_particle);

    Rep = denf * sqrt((uf - w_particle)^2 + v_particle^2) * dp / visf;
    if Rep < 800
        Cd = (24 / Rep) * (1 + (0.15 * Rep^(0.687)));
    else 
        Cd = 0.44;
    end

    Fdrag_x = (0.5 * pi * denf * dp^2 * Cd * sqrt((uf - w_particle)^2 + v_particle^2) * (uf - w_particle));
    Fdrag_y = (0.5 * pi * denf * dp^2 * Cd * sqrt((uf - w_particle)^2 + v_particle^2) * v_particle);

    F_tw = Fdrag_x / m_particle;
    F_tv = (Fg + Fb + Fdrag_y) / m_particle;

    dydt = [v_particle; F_tv; w_particle; F_tw];
end

% 4th order Runge-Kutta method for particle motion
function [y, x] = rungeKutta4(func, t, initial_conditions, yc, denf, u0, h, denp, dp, visf, g, Vp)
    dt = t(2) - t(1);
    y = zeros(length(t), 1);
    x = zeros(length(t), 1);
    
    y(1) = initial_conditions(1);
    v = initial_conditions(2);
    x(1) = initial_conditions(3);
    w = initial_conditions(4);

    for i = 1:numel(t)-1
        if y(i) < yc
            break; % Stop loop when y(i) < yc
        end

        k1 = func(t(i), [y(i), v, x(i), w], yc, denf, u0, h, denp, dp, visf, g, Vp);
        k2 = func(t(i) + dt/2, [y(i), v, x(i), w] + dt/2 * k1', yc, denf, u0, h, denp, dp, visf, g, Vp);
        k3 = func(t(i) + dt/2, [y(i), v, x(i), w] + dt/2 * k2', yc, denf, u0, h, denp, dp, visf, g, Vp);
        k4 = func(t(i) + dt, [y(i), v, x(i), w] + dt * k3', yc, denf, u0, h, denp, dp, visf, g, Vp);

        y(i+1) = y(i) + dt/6 * (k1(1) + 2*k2(1) + 2*k3(1) + k4(1));
        v = v + dt/6 * (k1(2) + 2*k2(2) + 2*k3(2) + k4(2));
        x(i+1) = x(i) + dt/6 * (k1(3) + 2*k2(3) + 2*k3(3) + k4(3));
        w = w + dt/6 * (k1(4) + 2*k2(4) + 2*k3(4) + k4(4));
    end

    y = y(1:i); % Particle y positions
    x = x(1:i); % Particle x positions
end

I have tried playing around with the code. I am only trying to learn this now suing a tutorial and this was one of the problems listed there with no solution. I am unfamiliar with most of the concepts and trying to learn so sometimes I miss out understanding. So even pointers to where I am missing stuff would be great.