import qualified Data.Stream as S
import Data.Stream (Stream(..))

-- Approximates c + integral f(t) dt for t=0 to infinity using simple Euler steps
-- Both the function and the steps are provided as streams
eulerIntegral :: Double -> Stream Double -> Stream Double -> Stream Double
eulerIntegral c ft dt = S.scan' (+) c $ S.zipWith (*) ft dt

-- Since scan' is strict it doesn't behave well regarding fixed points, 
-- e.g. the computation of the approximation of e blocks
e = es S.!! (10^6)
  where
    es = eulerIntegral 1 es dts
    dts = S.repeat 1e-6

-- Replacing S.scan' with myscan' fixes it
myscan' f b ~(a `Cons` as) = b `seq` (b `Cons` myscan' f (f b a) as)