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

Globals2D;

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

filter = Filter2D(0.5 * N,32);

rho = zeros(N_p,K);
u = zeros(N_p,K);
v = zeros(N_p,K);
p = zeros(N_p,K);

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

% 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 = 5;
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_sponge(rho,u,v,p,time_RK);
        
        % 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;
        
		rho = rho + rk4b(i) * rho_temp;
        u = u + rk4b(i) * u_temp;
        v = v + rk4b(i) * v_temp;
        p = p + rk4b(i) * 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