We introduce a new Neumann problem for the fractional Laplacian arising from a simple probabilistic consideration, and we discuss the basic properties of this model. We can consider both elliptic and parabolic equations in any domain. In addition, we formulate problems with nonhomogeneous Neumann conditions, and also with mixed Dirichlet and Neumann conditions, all of them having a clear probabilistic interpretation. We prove that solutions to the fractional heat equation with homogeneous Neumann conditions have the following natural properties: conservation of mass inside \Omega , decreasing energy, and convergence to a constant as t\to \infty . Moreover, for the elliptic case we give the variational formulation of the problem, and establish existence of solutions. We also study the limit properties and the boundary behavior induced by this nonlocal Neumann condition. For concreteness, one may think that our nonlocal analogue of the classical Neumann condition \partial_\nu u=0 on \partial\Omega consists in the nonlocal prescription \int_\Omega \frac{u(x)-u(y)}{|x-y|^{n+2s}}\,dy=0 \quad {\text{for }} x\in\mathbb R^n\setminus\overline{\Omega}. We made an effort to keep all the arguments at the simplest possible technical level, in order to clarify the connections between the different scientific fields that are naturally involved in the problem, and make the paper accessible also to a wide, non-specialistic public (for this scope, we also tried to use and compare different concepts and notations in a somehow more unified way).