IIRC it is not possible to express the indefinite integral in closed form. There is no doubt an infinite series that can be found. For certain values of the limits it is possible to evaluate the definite integral to an exact value.
It's not possible to antidifferentiate this function, but there is a trick you can use for a definite integral involving this function. If, for example, your integral is int(e^(x^2))dx, x = 0 to x = 1, this has the same value as int(e^(y^2))dy, y = 0 to y =1. If you multiply these together, the result will be the square of the number you want. The product of the two integrals is [int(e^(x^2))dx] * [int(e^(y^2))dy], which equals the double integral of e^(x^2 + y^2) dA over the region R = {(x, y) | 0 <= x <= 1, 0 <= y <= 1}. This double integral can be written as an iterated integral in polar form, as ∫∫e^(r^2)r dr dθ, where the limits of integration on the inner integral are 0 and the value of r that corresponds to x = 1 on the region R that was described earlier. (I'm leaving out some of the story here, since you also need the value of r along the upper edge of the rectangle R. Alternatively you can just integrate over the triangle with vertices (0,0), (1,0) and (1, 1) and double the result you get, using symmetry.) The limits of integration on the outer integral are θ = 0 and θ = π/2. Conversion to polar form introduces a factor of r when you go from dx dy to r dr dθ, which makes it possible to find an antiderivative for e^(r^2). Mark
I believe you can also evaluate this by finding the fourier transform and evaluating that at omega = 0.