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

Globals2D;

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

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
% zobecneni Advection2D Hesthaven - jeste stabilni
GLnodes = GaussLobattoNodes(N_fp); r_min = abs(GLnodes(2) - GLnodes(1));
dtscale = dtscale2D(); C = 1; dt = C * min(dtscale) * (2 / 3) * r_min / (a_0 + M_0); % Hesthaven volil C = 0.75
% jako Hesthaven v CFD
% dt = 1 / max(0.5 * ((N + 1) ^ 2) * Fscale(:) * (a_0 + M_0));
N_steps = ceil(T / dt); dt = T / N_steps;

current_time = 0;
for timestep = 1:N_steps
	for i = 1:5
		time_RK = current_time + rk4c(i) * dt;
		[rho_rhs,u_rhs,v_rhs,p_rhs] = ODEs_RHS(rho,u,v,p);
        
		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
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