// Time evolution
#pragma omp parallel default(none) shared(P, D, dt, dx, file, time, std::cout)
  {
    for (int step = 1; step <= numSteps; step++) {

    // Compute next time point
      diffusion1d_timestep(P, D, dt, dx);

    // Update time only on one thread
#pragma single
      {
        time += dt;

        // Periodically add data to the file
        if (step % outputEvery == 0 and step > 0)
          diffusion1d_output(file, step, time, P, outputcols);
      }
    }
  }