function [rho_rhs,u_rhs,v_rhs,w_rhs,p_rhs] = ODEs_RHS(rho,u,v,w,p)

Globals3D;

% local Lax-Friedrichs flux - strong form - fact that u_0 = v_0 = 0 used for simplification

% computing fluxes on the boundaries within the elements
rho_flux = u(vmapM) .* n_x(:) + v(vmapM) .* n_y(:) + w(vmapM) .* n_z(:);
u_flux = p(vmapM) .* n_x(:);
v_flux = p(vmapM) .* n_y(:);
w_flux = p(vmapM) .* n_z(:);
p_flux = u(vmapM) .* n_x(:) + v(vmapM) .* n_y(:) + w(vmapM) .* n_z(:);

% computing numerical fluxes
rho_numerical_flux = ((u(vmapM) + u(vmapP)) .* n_x(:) + (v(vmapM) + v(vmapP)) .* n_y(:) + (w(vmapM) + w(vmapP)) .* n_z(:) + (rho(vmapM) - rho(vmapP))) / 2;
u_numerical_flux = ((p(vmapM) + p(vmapP)) .* n_x(:) + (u(vmapM) - u(vmapP))) / 2;
v_numerical_flux = ((p(vmapM) + p(vmapP)) .* n_y(:) + (v(vmapM) - v(vmapP))) / 2;
w_numerical_flux = ((p(vmapM) + p(vmapP)) .* n_z(:) + (w(vmapM) - w(vmapP))) / 2;
p_numerical_flux = ((u(vmapM) + u(vmapP)) .* n_x(:) + (v(vmapM) + v(vmapP)) .* n_y(:) + (w(vmapM) + w(vmapP)) .* n_z(:) + (p(vmapM) - p(vmapP))) / 2;

% imposing boundary condition (unlimited domain) - postup Toulorge zalozeny na diagonalizaci matice A = n_x * A_1 + n_y * A_2
%for i = 1:length(mapB)
%   node_position_faces = mapB(i);
%   node_position_total = vmapB(i);
%   variables_vector = [rho(node_position_total); u(node_position_total); v(node_position_total); p(node_position_total)];
%   numerical_flux_face = TDplusTinv_flux_cell{node_position_faces} * variables_vector;
%   rho_numerical_flux(node_position_faces) = numerical_flux_face(1);
%   u_numerical_flux(node_position_faces) = numerical_flux_face(2);
%   v_numerical_flux(node_position_faces) = numerical_flux_face(3);
%   p_numerical_flux(node_position_faces) = numerical_flux_face(4);
%end
% imposing wall boundary condition - replacement for NRBC
scalar_product = u(vmapB) .* n_x(mapB) + v(vmapB) .* n_y(mapB) + w(vmapB) .* n_z(mapB);
u_modified = u(vmapB) - scalar_product .* n_x(mapB);
v_modified = v(vmapB) - scalar_product .* n_y(mapB);
w_modified = w(vmapB) - scalar_product .* n_z(mapB);
rho_numerical_flux(mapB) = u_modified .* n_x(mapB) + v_modified .* n_y(mapB) + w_modified .* n_z(mapB);
u_numerical_flux(mapB) = p(vmapB) .* n_x(mapB);
v_numerical_flux(mapB) = p(vmapB) .* n_y(mapB);
w_numerical_flux(mapB) = p(vmapB) .* n_z(mapB);
p_numerical_flux(mapB) = u_modified .* n_x(mapB) + v_modified .* n_y(mapB) + w_modified .* n_z(mapB);

% resizing both normal and numerical fluxes to appropiate arrays
rho_flux = reshape(rho_flux,N_faces * N_fp,K);
u_flux = reshape(u_flux,N_faces * N_fp,K);
v_flux = reshape(v_flux,N_faces * N_fp,K);
w_flux = reshape(w_flux,N_faces * N_fp,K);
p_flux = reshape(p_flux,N_faces * N_fp,K);

rho_numerical_flux = reshape(rho_numerical_flux,N_faces * N_fp,K);
u_numerical_flux = reshape(u_numerical_flux,N_faces * N_fp,K);
v_numerical_flux = reshape(v_numerical_flux,N_faces * N_fp,K);
w_numerical_flux = reshape(w_numerical_flux,N_faces * N_fp,K);
p_numerical_flux = reshape(p_numerical_flux,N_faces * N_fp,K);

% building the right-hand side of the ODE system
% [u_x,u_y,u_z] = Grad3D(u); [v_x,v_y,v_z] = Grad3D(v); [w_x,w_y,w_z] = Grad3D(v); [p_x,p_y,p_z] = Grad3D(p);
u_r = D_r * u; u_s = D_s * u; u_t = D_t * u; u_x = r_x .* u_r + s_x .* u_s + t_x .* u_t;
v_r = D_r * v; v_s = D_s * v; v_t = D_t * v; v_y = r_y .* v_r + s_y .* v_s + t_y .* v_t;
w_r = D_r * w; w_s = D_s * w; w_t = D_t * w; w_z = r_z .* w_r + s_z .* w_s + t_z .* w_t;
p_r = D_r * p; p_s = D_s * p; p_t = D_t * p; p_x = r_x .* p_r + s_x .* p_s + t_x .* p_t; p_y = r_y .* p_r + s_y .* p_s + t_y .* p_t; p_z = r_z .* p_r + s_z .* p_s + t_z .* p_t;
rho_rhs = -u_x - v_y - w_z + LIFT * (Fscale .* (rho_flux - rho_numerical_flux));
u_rhs = -p_x + LIFT * (Fscale .* (u_flux - u_numerical_flux));
v_rhs = -p_y + LIFT * (Fscale .* (v_flux - v_numerical_flux));
w_rhs = -p_z + LIFT * (Fscale .* (w_flux - w_numerical_flux));
p_rhs = -u_x - v_y - w_z  + LIFT * (Fscale .* (p_flux - p_numerical_flux));

end