function u_rhs = ODEs_RHS(u,time_RK)

Globals3D;

% local Lax-Friedrichs flux - strong form

% computing flux on the boundaries within the elements
flux = (a_advec * u(vmapM)) .* n_x(:) + (b_advec * u(vmapM)) .* n_y(:) + (c_advec * u(vmapM)) .* n_z(:);

% computing numerical flux
numerical_flux = ((a_advec * (u(vmapM) + u(vmapP))) .* n_x(:) + (b_advec * (u(vmapM) + u(vmapP))) .* n_y(:) + (c_advec * (u(vmapM) + u(vmapP))) .* n_z(:) + (u(vmapM) - u(vmapP)) .* abs(a_advec * n_x(:) + b_advec * n_y(:) + c_advec * n_z(:))) / 2;

% imposing inflow condition
x_in = x(vmapI); y_in = y(vmapI); z_in = z(vmapI); u_inflow = BC(x_in,y_in,z_in,time_RK);
numerical_flux(mapI) = ((a_advec * (u(vmapI) + u_inflow)) .* n_x(mapI) + (b_advec * (u(vmapI) + u_inflow)) .* n_y(mapI) + (c_advec * (u(vmapI) + u_inflow)) .* n_z(mapI) + (u(vmapI) - u_inflow) .* abs(a_advec * n_x(mapI) + b_advec * n_y(mapI) + c_advec * n_z(mapI))) / 2;

% imposing outflow condition - we extrapolate u - therefore the numerical flux takes form
numerical_flux(mapO) = (a_advec * u(vmapO)) .* n_x(mapO) + (b_advec * u(vmapO)) .* n_y(mapO) + (c_advec * u(vmapO)) .* n_z(mapO);

% resizing both normal and numerical fluxes to appropiate arrays
flux = reshape(flux,N_fp * N_faces,K);
numerical_flux = reshape(numerical_flux,N_fp * N_faces,K);

% building the right-hand side of the ODE system
[u_x,u_y,u_z] = Grad3D(u);
u_rhs = -a_advec * u_x - b_advec * u_y - c_advec * u_z + LIFT * (Fscale .* (flux - numerical_flux));

end