function [rho_rhs,u_rhs,v_rhs,p_rhs] = ODEs_RHS_PML(rho,u,v,p,q_rho,q_u,q_v,q_p)

Globals2D;

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

% computing fluxes on the boundaries within the elements
rho_flux = (u_0 * rho(vmapM) + rho_0 * u(vmapM)) .* n_x(:) + (rho_0 * v(vmapM)) .* n_y(:);
u_flux = (u_0 * u(vmapM) + (1 / rho_0) * p(vmapM)) .* n_x(:);
v_flux = (u_0 * v(vmapM)) .* n_x(:) + ((1 / rho_0) * p(vmapM)) .* n_y(:);
p_flux = (poisson_const * p_0 * u(vmapM) + u_0 * p(vmapM)) .* n_x(:) + (poisson_const * p_0 * v(vmapM)) .* n_y(:);

rho_flux(mapPML_M) = rho_flux(mapPML_M) + (rho_0 * q_v(vmapPML_M)) .* sigma_x_PML_times_n_y(mapPML_M);
% no term to be added in the second equation
v_flux(mapPML_M) = v_flux(mapPML_M) + ((1 / rho_0) * q_p(vmapPML_M)) .* sigma_x_PML_times_n_y(mapPML_M);
p_flux(mapPML_M) = p_flux(mapPML_M) + (poisson_const * p_0 * q_v(vmapPML_M)) .* sigma_x_PML_times_n_y(mapPML_M);

% computing numerical fluxes
rho_numerical_flux = ((u_0 * (rho(vmapM) + rho(vmapP)) + rho_0 * (u(vmapM) + u(vmapP))) .* n_x(:) + (rho_0 * (v(vmapM) + v(vmapP))) .* n_y(:) + (rho(vmapM) - rho(vmapP)) .* ((1 - alpha_flux) * C_flux(:))) / 2;
u_numerical_flux = ((u_0 * (u(vmapM) + u(vmapP)) + (1 / rho_0) * (p(vmapM) + p(vmapP))) .* n_x(:) + (u(vmapM) - u(vmapP)) .* ((1 - alpha_flux) * C_flux(:))) / 2;
v_numerical_flux = ((u_0 * (v(vmapM) + v(vmapP))) .* n_x(:) + ((1 / rho_0) * (p(vmapM) + p(vmapP))) .* n_y(:) + (v(vmapM) - v(vmapP)) .* ((1 - alpha_flux) * C_flux(:))) / 2;
p_numerical_flux = ((poisson_const * p_0 * (u(vmapM) + u(vmapP)) + u_0 * (p(vmapM) + p(vmapP))) .* n_x(:) + (poisson_const * p_0 * (v(vmapM) + v(vmapP))) .* n_y(:) + (p(vmapM) - p(vmapP)) .* ((1 - alpha_flux) * C_flux(:))) / 2;

rho_numerical_flux(mapPML_M) = rho_numerical_flux(mapPML_M) + (rho_0 * (q_v(vmapPML_M) + q_v(vmapPML_P))) .* sigma_x_PML_times_n_y(mapPML_M) / 2;
% no term to be added in the second equation
v_numerical_flux(mapPML_M) = v_numerical_flux(mapPML_M) + ((1 / rho_0) * (q_p(vmapPML_P) + q_p(vmapPML_P))) .* sigma_x_PML_times_n_y(mapPML_M) / 2;
p_numerical_flux(mapPML_M) = p_numerical_flux(mapPML_M) + (poisson_const * p_0 * (q_v(vmapPML_M) + q_v(vmapPML_P))) .* sigma_x_PML_times_n_y(mapPML_M) / 2;

% imposing boundary condition (solid wall) - neslo by toto napsat vektorove, tj. efektivneji nez pres cyklus?
for i = 1:length(mapW)
    node_position_faces = mapW(i);
    node_position_total = vmapW(i);
    velocities_modified = [u(node_position_total); v(node_position_total)] - (u(node_position_total) * n_x(node_position_faces) + v(node_position_total) * n_y(node_position_faces)) * [n_x(node_position_faces); n_y(node_position_faces)];
    variables_vector = [rho(node_position_total); velocities_modified; p(node_position_total)];
    numerical_flux_face = (n_x(node_position_faces) * A_1 + n_y(node_position_faces) * A_2) * 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
% rho_numerical_flux(mapW_PML) = rho_numerical_flux(mapW_PML) + (sigma_x_PML(vmapW_PML) .* (rho_0 * q_v(vmapW_PML))) .* n_y(mapW_PML);
% % no term to be added in the second equation
% v_numerical_flux(mapW_PML) = v_numerical_flux(mapW_PML) + (sigma_x_PML(vmapW_PML) .* ((1 / rho_0) * q_p(vmapW_PML))) .* n_y(mapW_PML);
% p_numerical_flux(mapW_PML) = p_numerical_flux(mapW_PML) + (sigma_x_PML(vmapW_PML) .* (poisson_const * p_0 * q_v(vmapW_PML))) .* n_y(mapW_PML);

% closing the PML with Dirichlet BC, viz. Hu 2001
% inflow
nonPMLpart = (rho_0 * u(vmapI)) .* n_x(mapI);
PMLpart = (rho_0 * q_v(vmapI)) .* sigma_x_PML_times_n_y(mapI);
rho_numerical_flux(mapI) = nonPMLpart + PMLpart;
nonPMLpart = (u_0 * u(vmapI)) .* n_x(mapI);
%PMLpart = zeros(length(mapI),1);
%u_numerical_flux(mapI) = nonPMLpart + PMLpart; % no PML part
u_numerical_flux(mapI) = nonPMLpart;
%nonPMLpart = zeros(length(mapI),1);
PMLpart = ((1 / rho_0) * q_p(vmapI)) .* sigma_x_PML_times_n_y(mapI);
%v_numerical_flux(mapI) = nonPMLpart + PMLpart; % no nonPML part
v_numerical_flux(mapI) = PMLpart;
nonPMLpart = (poisson_const * p_0 * u(vmapI)) .* n_x(mapI);
PMLpart = (poisson_const * p_0 * q_v(vmapI)) .* sigma_x_PML_times_n_y(mapI);
p_numerical_flux(mapI) = nonPMLpart + PMLpart;
% outflow
nonPMLpart = (u_0 * rho(vmapO) + rho_0 * u(vmapO)) .* n_x(mapO) + (rho_0 * v(vmapO)) .* n_y(mapO);
PMLpart = (rho_0 * q_v(vmapO)) .* sigma_x_PML_times_n_y(mapO);
rho_numerical_flux(mapO) = nonPMLpart + PMLpart;
nonPMLpart = (u_0 * u(vmapO)) .* n_x(mapO);
%PMLpart = zeros(length(mapO),1);
%u_numerical_flux(mapO) = nonPMLpart + PMLpart; % no PML part
u_numerical_flux(mapO) = nonPMLpart;
nonPMLpart = (u_0 * v(vmapO)) .* n_x(mapO);
PMLpart = ((1 / rho_0) * q_p(vmapO)) .* sigma_x_PML_times_n_y(mapO);
v_numerical_flux(mapO) = nonPMLpart + PMLpart;
nonPMLpart = (poisson_const * p_0 * u(vmapO)) .* n_x(mapO) + (poisson_const * p_0 * v(vmapO)) .* n_y(mapO);
PMLpart = (poisson_const * p_0 * q_v(vmapO)) .* sigma_x_PML_times_n_y(mapO);
p_numerical_flux(mapO) = nonPMLpart + PMLpart;

% 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);
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);
p_numerical_flux = reshape(p_numerical_flux,N_faces * N_fp,K);

% building the right-hand side of the ODE system
[rho_x,rho_y] = Grad2D(rho); [u_x,u_y] = Grad2D(u); [v_x,v_y] = Grad2D(v); [p_x,p_y] = Grad2D(p);
rho_rhs = -u_0 * rho_x - rho_0 * (u_x + v_y) + LIFT * (Fscale .* (rho_flux - rho_numerical_flux));
u_rhs = -u_0 * u_x - (1 / rho_0) * p_x + LIFT * (Fscale .* (u_flux - u_numerical_flux));
v_rhs = -u_0 * v_x - (1 / rho_0) * p_y + LIFT * (Fscale .* (v_flux - v_numerical_flux));
p_rhs = -poisson_const * p_0 * (u_x + v_y) - u_0 * p_x + LIFT * (Fscale .* (p_flux - p_numerical_flux));

[q_v_x,q_v_y] = Grad2D(q_v); [q_p_x,q_p_y] = Grad2D(q_p); % q_rho a q_u neni potreba derivovat
derivatives = sigma_x_PML(:,elements_PML) .* (rho_0 * q_v_y(:,elements_PML));
non_derivatives = sigma_x_PML(:,elements_PML) .* rho(:,elements_PML);
stabilization = coef_PML * (sigma_x_PML(:,elements_PML) .* (u_0 * rho(:,elements_PML) + rho_0 * u(:,elements_PML)));
rho_rhs(:,elements_PML) = rho_rhs(:,elements_PML) - derivatives - non_derivatives - stabilization; 
%derivatives = zeros(N_p,length(elements_PML));
non_derivatives = sigma_x_PML(:,elements_PML) .* u(:,elements_PML);
stabilization = coef_PML * (sigma_x_PML(:,elements_PML) .* (u_0 * u(:,elements_PML) + (1 / rho_0) * p(:,elements_PML)));
%u_rhs(:,elements_PML) = u_rhs(:,elements_PML) - derivatives - non_derivatives - stabilization; % no derivatives part
u_rhs(:,elements_PML) = u_rhs(:,elements_PML) - non_derivatives - stabilization;
derivatives = sigma_x_PML(:,elements_PML) .* ((1 / rho_0) * q_p_y(:,elements_PML));
non_derivatives = sigma_x_PML(:,elements_PML) .* v(:,elements_PML);
stabilization = coef_PML * (sigma_x_PML(:,elements_PML) .* (u_0 * v(:,elements_PML)));
v_rhs(:,elements_PML) = v_rhs(:,elements_PML) - derivatives - non_derivatives - stabilization;
derivatives = sigma_x_PML(:,elements_PML) .* (poisson_const * p_0 * q_v_y(:,elements_PML));
non_derivatives = sigma_x_PML(:,elements_PML) .* p(:,elements_PML);
stabilization = coef_PML * (sigma_x_PML(:,elements_PML) .* (poisson_const * p_0 * u(:,elements_PML) + u_0 * p(:,elements_PML)));
p_rhs(:,elements_PML) = p_rhs(:,elements_PML) - derivatives - non_derivatives - stabilization;

end