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 = ['Bajko_Problem3_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;
%N_steps = T / dt; dt = T / N_steps;
interval = 10;
N_steps_states = ceil(interval / dt); dt = interval / N_steps_states;
N_steps = ceil(T / dt);


% reseni na hranici
% (x,y) = (100,0)
x_boundary_R = 100; y_boundary_R = 0;
rho_boundary_R = zeros(1,N_steps + 1);
rho_boundary_R(1) = rho_IC(x_boundary_R,y_boundary_R);
u_boundary_R = zeros(1,N_steps + 1);
u_boundary_R(1) = u_IC(x_boundary_R,y_boundary_R);
v_boundary_R = zeros(1,N_steps + 1);
v_boundary_R(1) = v_IC(x_boundary_R,y_boundary_R);
p_boundary_R = zeros(1,N_steps + 1);
p_boundary_R(1) = p_IC(x_boundary_R,y_boundary_R);
for element_boundary_R = 1:K
    vertices = P(:,E(:,element_boundary_R));
    A = [vertices(1,2) - vertices(1,1), vertices(1,3) - vertices(1,1); vertices(2,2) - vertices(2,1), vertices(2,3) - vertices(2,1)];
    b = [x_boundary_R - vertices(1,1); y_boundary_R - vertices(2,1)];
    koef_temp = A\b;
    koef = [1 - sum(koef_temp); koef_temp];
    
    if (min(koef) >= 0)
        break;
    end
end
T_matrix = [vertices(1,1) - vertices(1,3), vertices(1,2) - vertices(1,3); vertices(2,1) - vertices(2,3), vertices(2,2) - vertices(2,3)];
det_T = det(T_matrix);
lambda2 = ((vertices(2,2) - vertices(2,3)) * (x_boundary_R - vertices(1,3)) + (vertices(1,3) - vertices(1,2)) * (y_boundary_R - vertices(2,3))) / det_T;
lambda3 = ((vertices(2,3) - vertices(2,1)) * (x_boundary_R - vertices(1,3)) + (vertices(1,1) - vertices(1,3)) * (y_boundary_R - vertices(2,3))) / det_T;
lambda1 = 1 - lambda2 - lambda3;
r = -lambda2 + lambda3 - lambda1;
s = -lambda2 - lambda3 + lambda1;
[a,b] = rstoab(r,s);
values_mode = zeros(N_p,1);
temp = 1;
for k = 0:N
    for l = 0:(N - k)
        values_mode(temp,:) = ModalBaseValues2D(a,b,k,l);
        temp = temp + 1;
    end
end
values_node = inv(V') * values_mode;
% (x,y) = (-100,0)
x_boundary_L = -100; y_boundary_L = 0;
rho_boundary_L = zeros(1,N_steps + 1);
rho_boundary_L(1) = rho_IC(x_boundary_L,y_boundary_L);
u_boundary_L = zeros(1,N_steps + 1);
u_boundary_L(1) = u_IC(x_boundary_L,y_boundary_L);
v_boundary_L = zeros(1,N_steps + 1);
v_boundary_L(1) = v_IC(x_boundary_L,y_boundary_L);
p_boundary_L = zeros(1,N_steps + 1);
p_boundary_L(1) = p_IC(x_boundary_L,y_boundary_L);
for element_boundary_L = 1:K
    vertices = P(:,E(:,element_boundary_L));
    A = [vertices(1,2) - vertices(1,1), vertices(1,3) - vertices(1,1); vertices(2,2) - vertices(2,1), vertices(2,3) - vertices(2,1)];
    b = [x_boundary_L - vertices(1,1); y_boundary_L - vertices(2,1)];
    koef_temp = A\b;
    koef = [1 - sum(koef_temp); koef_temp];
    
    if (min(koef) >= 0)
        break;
    end
end
T_matrix = [vertices(1,1) - vertices(1,3), vertices(1,2) - vertices(1,3); vertices(2,1) - vertices(2,3), vertices(2,2) - vertices(2,3)];
det_T = det(T_matrix);
lambda2 = ((vertices(2,2) - vertices(2,3)) * (x_boundary_L - vertices(1,3)) + (vertices(1,3) - vertices(1,2)) * (y_boundary_L - vertices(2,3))) / det_T;
lambda3 = ((vertices(2,3) - vertices(2,1)) * (x_boundary_L - vertices(1,3)) + (vertices(1,1) - vertices(1,3)) * (y_boundary_L - vertices(2,3))) / det_T;
lambda1 = 1 - lambda2 - lambda3;
r = -lambda2 + lambda3 - lambda1;
s = -lambda2 - lambda3 + lambda1;
[a,b] = rstoab(r,s);
values_mode = zeros(N_p,1);
temp = 1;
for k = 0:N
    for l = 0:(N - k)
        values_mode(temp,:) = ModalBaseValues2D(a,b,k,l);
        temp = temp + 1;
    end
end
values_node = inv(V') * values_mode;


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

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 (timestep ~= 0) && (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
    
    % solution on boundary
    % (x,y) = (100,0)
    coef = rho(:,element_boundary_R)';
    rho_boundary_R(timestep + 1) = coef * values_node;
    coef = u(:,element_boundary_R)';
    u_boundary_R(timestep + 1) = coef * values_node;
    coef = v(:,element_boundary_R)';
    v_boundary_R(timestep + 1) = coef * values_node;
    coef = p(:,element_boundary_R)';
    p_boundary_R(timestep + 1) = coef * values_node;
    % (x,y) = (-100,0)
    coef = rho(:,element_boundary_L)';
    rho_boundary_L(timestep + 1) = coef * values_node;
    coef = u(:,element_boundary_L)';
    u_boundary_L(timestep + 1) = coef * values_node;
    coef = v(:,element_boundary_L)';
    v_boundary_L(timestep + 1) = coef * values_node;
    coef = p(:,element_boundary_L)';
    p_boundary_L(timestep + 1) = coef * values_node;
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