function [rho,u,v,p] = LEETimestepping2D_Problem5_PML(rho,u,v,p,T)

Globals2D;

mesh_size_string = sprintf('%0.2f',mesh_size_parametr);
mesh_size_string = strrep(mesh_size_string,'.',',');
filename = ['PML_test_Results_in_time_for_h=',mesh_size_string,',_N=',num2str(N),'.txt'];
fileID = fopen(filename,'w');

filter = CutOffFilter2D(N,0.95);

rho_temp = zeros(N_p,K);
u_temp = zeros(N_p,K);
v_temp = zeros(N_p,K);
p_temp = zeros(N_p,K);

% arrays for auxiliary variables in q
q_rho = zeros(N_p,K);
q_u = zeros(N_p,K);
q_v = zeros(N_p,K);
q_p = zeros(N_p,K);

q_rho_temp = zeros(N_p,K);
q_u_temp = zeros(N_p,K);
q_v_temp = zeros(N_p,K);
q_p_temp = zeros(N_p,K);

% coef for one term in PML formulation preprocessed
coef_PML = u_0 / (1 - u_0 ^ 2);

% time-step choice for stable scheme
%GLnodes = GaussLobattoNodes(N_fp); r_min = abs(GLnodes(2) - GLnodes(1));
%dtscale = dtscale2D(); C = 0.75; dt = C * min(dtscale) * (2 / 3) * r_min / (a_0 + M_0);
dt = 1 / max(0.5 * ((N + 1) ^ 2) * Fscale(:) * (a_0 + M_0)); dt = 0.9 * dt;
%N_steps = ceil(T / dt); dt = T / N_steps;
%N_steps = T / dt; dt = T / N_steps;
interval = 1;
N_steps_states = ceil(interval / dt); dt = interval / N_steps_states;
N_steps = ceil(T / dt);

states_counter = 1;
p_states = cell(1,T / interval);

current_time = 0;
for timestep = 1:N_steps
	for i = 1:5 % 5 stages of low-storage Runge-Kutta method
		time_RK = current_time + rk4c(i) * dt;
        [rho_rhs,u_rhs,v_rhs,p_rhs] = ODEs_RHS_PML(rho,u,v,p,q_rho,q_u,q_v,q_p);
        
        % filtering the computed right-hand sides of the ODE system
        rho_rhs = filter * rho_rhs;
        u_rhs = filter * u_rhs;
        v_rhs = filter * v_rhs;
        p_rhs = filter * p_rhs;
        
		rho_temp = rk4a(i) * rho_temp + dt * rho_rhs;
        u_temp = rk4a(i) * u_temp + dt * u_rhs;
        v_temp = rk4a(i) * v_temp + dt * v_rhs;
        p_temp = rk4a(i) * p_temp + dt * p_rhs;
        q_rho_temp = rk4a(i) * q_rho_temp + dt * rho; % q_rho_rhs = rho
        q_u_temp = rk4a(i) * q_u_temp + dt * u; % q_u_rhs = u
        q_v_temp = rk4a(i) * q_v_temp + dt * v; % q_v_rhs = v
        q_p_temp = rk4a(i) * q_p_temp + dt * p; % q_p_rhs = p
        
		rho = rho + rk4b(i) * rho_temp;
        u = u + rk4b(i) * u_temp;
        v = v + rk4b(i) * v_temp;
        p = p + rk4b(i) * p_temp;
        q_rho = q_rho + rk4b(i) * q_rho_temp;
        q_u = q_u + rk4b(i) * q_u_temp;
        q_v = q_v + rk4b(i) * q_v_temp;
        q_p = q_p + rk4b(i) * q_p_temp;
	end
	current_time = current_time + dt;
    
    if (rem(timestep,100) == 0)
        fprintf(2,'%i of total %i timesteps processed.\n',timestep,N_steps);
        
        fprintf(fileID,'Time %0.3fs:\n',current_time);
        fprintf(fileID,'Min rho: %e, \t max rho: %e\n',min(min(rho)),max(max(rho)));
        fprintf(fileID,'Min u: %e, \t max u: %e\n',min(min(u)),max(max(u)));
        fprintf(fileID,'Min v: %e, \t max v: %e\n',min(min(v)),max(max(v)));
        fprintf(fileID,'Min p: %e, \t max p: %e\n\n',min(min(p)),max(max(p)));
    end
    
    if (rem(timestep,N_steps_states) == 0)
        fprintf(2,'Min p: %e, \t max p: %e\n',min(min(p)),max(max(p)));
        
        p_states{states_counter} = p;
        states_counter = states_counter + 1;
    end
end

fprintf(fileID,'Time %0.3fs:\n',current_time);
fprintf(fileID,'Min rho: %e, \t max rho: %e\n',min(min(rho)),max(max(rho)));
fprintf(fileID,'Min u: %e, \t max u: %e\n',min(min(u)),max(max(u)));
fprintf(fileID,'Min v: %e, \t max v: %e\n',min(min(v)),max(max(v)));
fprintf(fileID,'Min p: %e, \t max p: %e\n\n',min(min(p)),max(max(p)));
        
fclose(fileID);

end