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 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(:);
u_flux = p(vmapM) .* n_x(:);
v_flux = p(vmapM) .* n_y(:);
p_flux = u(vmapM) .* n_x(:) + v(vmapM) .* n_y(:);

rho_flux(mapPML_M) = rho_flux(mapPML_M) + q_u(vmapPML_M) .* sigma_y_PML_times_n_x(mapPML_M) + q_v(vmapPML_M) .* sigma_x_PML_times_n_y(mapPML_M);
u_flux(mapPML_M) = u_flux(mapPML_M) + q_p(vmapPML_M) .* sigma_y_PML_times_n_x(mapPML_M);
v_flux(mapPML_M) = v_flux(mapPML_M) + q_p(vmapPML_M) .* sigma_x_PML_times_n_y(mapPML_M);
p_flux(mapPML_M) = p_flux(mapPML_M) + q_u(vmapPML_M) .* sigma_y_PML_times_n_x(mapPML_M) + q_v(vmapPML_M) .* sigma_x_PML_times_n_y(mapPML_M);

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

rho_numerical_flux(mapPML_M) = rho_numerical_flux(mapPML_M) + ((q_u(vmapPML_M) + q_u(vmapPML_P)) .* sigma_y_PML_times_n_x(mapPML_M) + (q_v(vmapPML_M) + q_v(vmapPML_P)) .* sigma_x_PML_times_n_y(mapPML_M)) / 2;
u_numerical_flux(mapPML_M) = u_numerical_flux(mapPML_M) + ((q_p(vmapPML_M) + q_p(vmapPML_P)) .* sigma_y_PML_times_n_x(mapPML_M)) / 2;
v_numerical_flux(mapPML_M) = v_numerical_flux(mapPML_M) + ((q_p(vmapPML_M) + q_p(vmapPML_P)) .* sigma_x_PML_times_n_y(mapPML_M)) / 2;
p_numerical_flux(mapPML_M) = p_numerical_flux(mapPML_M) + ((q_u(vmapPML_M) + q_u(vmapPML_P)) .* sigma_y_PML_times_n_x(mapPML_M) + (q_v(vmapPML_M) + q_v(vmapPML_P)) .* sigma_x_PML_times_n_y(mapPML_M)) / 2;

% closing the PML with characteristic-based BC - pro tuto implementaci je potreba dopocitat vlastni vektory
%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

% closing the PML with Dirichlet BC, viz. Hu 2001 - pouzita specifikace numerickeho toku pouhym dosazenim okr. podminek, tj. nebyl pouzit Lax-Friedrichs tok
% outflow - cela hranice je v tomto pripade outflow pro u_0 = v_0 = 0;
nonPMLpart = u(vmapB) .* n_x(mapB) + v(vmapB) .* n_y(mapB);
PMLpart = q_u(vmapB) .* sigma_y_PML_times_n_x(mapB) + q_v(vmapB) .* sigma_x_PML_times_n_y(mapB);
rho_numerical_flux(mapB) = nonPMLpart + PMLpart;
%nonPMLpart = zeros(length(mapB),1);
PMLpart = q_p(vmapB) .* sigma_y_PML_times_n_x(mapB);
%u_numerical_flux(mapB) = nonPMLpart + PMLpart; % no nonPML part
u_numerical_flux(mapB) = PMLpart;
%nonPMLpart = zeros(length(mapB),1);
PMLpart = q_p(vmapB) .* sigma_x_PML_times_n_y(mapB);
%v_numerical_flux(mapB) = nonPMLpart + PMLpart; % no nonPML part
v_numerical_flux(mapB) = PMLpart;
nonPMLpart = u(vmapB) .* n_x(mapB) + v(vmapB) .* n_y(mapB);
PMLpart = q_u(vmapB) .* sigma_y_PML_times_n_x(mapB) + q_v(vmapB) .* sigma_x_PML_times_n_y(mapB);
p_numerical_flux(mapB) = 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
[u_x,u_y] = Grad2D(u); [v_x,v_y] = Grad2D(v); [p_x,p_y] = Grad2D(p);
rho_rhs = -u_x - v_y + 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));
p_rhs = -u_x - v_y + LIFT * (Fscale .* (p_flux - p_numerical_flux));

[q_u_x,q_u_y] = Grad2D(q_u); [q_v_x,q_v_y] = Grad2D(q_v); [q_p_x,q_p_y] = Grad2D(q_p);
derivatives = sigma_y_PML(:,elements_PML) .* q_u_x(:,elements_PML) + sigma_x_PML(:,elements_PML) .* q_v_y(:,elements_PML);
non_derivatives = sigma_x_PML_plus_sigma_y_PML(:,elements_PML) .* rho(:,elements_PML) + sigma_x_PML_times_sigma_y_PML(:,elements_PML) .* q_rho(:,elements_PML);
rho_rhs(:,elements_PML) = rho_rhs(:,elements_PML) - derivatives - non_derivatives; 
derivatives = sigma_y_PML(:,elements_PML) .* q_p_x(:,elements_PML);
non_derivatives = sigma_x_PML_plus_sigma_y_PML(:,elements_PML) .* u(:,elements_PML) + sigma_x_PML_times_sigma_y_PML(:,elements_PML) .* q_u(:,elements_PML);
u_rhs(:,elements_PML) = u_rhs(:,elements_PML) - derivatives - non_derivatives;
derivatives = sigma_x_PML(:,elements_PML) .* q_p_y(:,elements_PML);
non_derivatives = sigma_x_PML_plus_sigma_y_PML(:,elements_PML) .* v(:,elements_PML) + sigma_x_PML_times_sigma_y_PML(:,elements_PML) .* q_v(:,elements_PML);
v_rhs(:,elements_PML) = v_rhs(:,elements_PML) - derivatives - non_derivatives;
derivatives = sigma_y_PML(:,elements_PML) .* q_u_x(:,elements_PML) + sigma_x_PML(:,elements_PML) .* q_v_y(:,elements_PML);
non_derivatives = sigma_x_PML_plus_sigma_y_PML(:,elements_PML) .* p(:,elements_PML) + sigma_x_PML_times_sigma_y_PML(:,elements_PML) .* q_p(:,elements_PML);
p_rhs(:,elements_PML) = p_rhs(:,elements_PML) - derivatives - non_derivatives;

end