function u = Timestepping(u,T)

Globals3D;

u_temp = zeros(N_p,K);

% time-step choice for stable scheme
GLnodes = GaussLobattoNodes(N + 1); r_min = abs(GLnodes(2) - GLnodes(1));
dtscale = dtscale3D(); C = 0.125; dt = C * min(dtscale) * (2 / 3) * r_min / sqrt(a_advec ^ 2 + b_advec ^ 2 + c_advec ^ 2);
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;
		u_rhs = ODEs_RHS(u,time_RK);
		u_temp = rk4a(i) * u_temp + dt * u_rhs;
		u = u + rk4b(i) * u_temp;
	end
	current_time = current_time + dt;
    
    if (rem(timestep,100) == 0)
        fprintf(2,'%i of total %i timesteps processed.\n',timestep,N_steps);
    end
end

end