An initial value problem for a first order differential equation with an unbounded constant operator coefficient $A$ in Hilbert space is considered. We give the definition of a $\sigma$-solution and using the Cayley transform we deduce an explicit formula for the solution in case the operator $A$ is self-adjoint and positively definite. On the basis of this formula we propose a numerical algorithm for the approximate solution of the initial value problem and give an error estimate. It turns out that, contrary to the case of a bounded operator $A$, the rate of convergence is not exponential but only polynomial and depends on the smoothness of the initial data. It is proved that the approximate solution is a best approximation in some Hilbert subspace. An example concerning the homogeneous heat equation is given.