Mike gives a reasonable approach. We could use the nth order approximation from sums.
Another would be to do a Taylor series approximation of the function and integrate this approximation from 0 to 1.
As the Taylor series for e^x is ∑x^n/n! , we substitute -x^2/5 in this expression.
Then, we have ∑(-1)^n x^(2n)/(5^n n!)
As we integrate this from 0 to 1, we get ∑(-1)^n /((2n+1)5^n n!) where the sum is for n = 0 to infinity.
Note that the absolute value of the ratio between successive terms is less than 1/(5n)
The term for n = 11 is -2.2 * 10^-17, which means that, as this is an alternating series, our error is less than this for n = 10.
Of course, we add 4 to the result.
It turns out that the cumulative distribution for the normal is on most calculators. The density function shown is that of the normal with standard deviation sqrt(5/2) with a necessary multiplication correction (as the normal density has a constant 1/(sqrt(2π)σ) of sqrt(5π)
Then, sqrt(5π) (norm.dist(1,0,sqrt(5/2),1)-.5) = 0.937150028797979, and 4 + 0.937150028797979 = 4.937150028797979
Our result at n = 10 matches this.