Numerically and analytically solve the radioactive decay problem (dN/dt = -kN) and plot the results.

euler_decay.hs

module Main where

import Graphics.EasyPlot

dt  = 0.001  -- Timestep
tau = 0.05   -- Decay constant
n0  = 1000   -- Initial number of particles

main = do
    plot X11 $ [Data2D [Title "Euler", Color Red, Style Lines] [] numerical,
                Data2D [Title "Actual", Color Blue, Style Lines] [] analytic]

-- Analytic results
analytic :: [(Double, Double)]
analytic = [(t, n0 * (exp (-t/tau))) | t <- [0.0, dt .. tau*5]]

-- Numerical results
numerical :: [(Double, Double)]
numerical = [(t, num t) | t <- [0.0, dt..tau*5]]

-- Numerically solve dN/dt = -kN
num :: Double -> Double
num t
    | t > 0     = let t' = num (t - dt) in t' - dt / tau * t'
    | otherwise = n0