function u_rhs = ODEs_RHS(u,time_RK)

Globals2D;

% 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(:);

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

% imposing inflow condition
u_inflow = BC(x(vmapI),y(vmapI),time_RK);
numerical_flux(mapI) = ((a_advec * (u(vmapI) + u_inflow)) .* n_x(mapI) + (b_advec * (u(vmapI) + u_inflow)) .* n_y(mapI) + (u(vmapI) - u_inflow) .* ((1 - alpha_flux) * abs(a_advec * n_x(mapI)+ b_advec * n_y(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);

% 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] = Grad2D(u);
u_rhs = -a_advec * u_x - b_advec * u_y + LIFT * (Fscale .* (flux - numerical_flux));

% verze z knihy - local Lax-Friedrichs flux
% du = (u(vmapM) - u(vmapP)) .* (a_advec * n_x(:) + b_advec * n_y(:) - (1 - alpha_flux) * abs(a_advec * n_x(:) + b_advec * n_y(:))) / 2;
% x_in = x(vmapI);
% y_in = y(vmapI);
% uin = BC(x_in,y_in,time_RK);
% du(mapI) = (u(vmapI) - uin) .* (a_advec * n_x(mapI) + b_advec * n_y(mapI) - (1 - alpha_flux) * abs(a_advec * n_x(mapI) + b_advec * n_y(mapI))) / 2;
% du(mapO) = 0;
% du = reshape(du,N_faces * N_fp,K);
% 
% [u_x,u_y] = Grad2D(u);
% u_rhs = -a_advec * u_x - b_advec * u_y + LIFT * (Fscale .* du);

end