To test accuracy of discrete approximation in time the equation of diffusion has been examined. To find **P(x,y,z,t)** particles distribution at t_{n} moment in time the Taylor expansion with β = 0 (the Euler ,,forward'' scheme) has been applied. The diffusion equation with the following initial condition g(., 0) = -x(x - π)y(y - π)z(z - π) for the cubic domain and g(., 0) = |(r - R_{0})z(z - π)| where R_{0} denotes a radius of the cylindrical domain has been solved. The boundary value of the **P(x,y,z,t)** is set as 0. Results for both domains: cubic (with uniform elements - see Fig.3a) and cylindrical (see Fig.1), this one tuned to the designed element volume set as 0.5 with the Metropolis recipe, are shown in Figs 5c, 6a and 6b. The diffusion coefficient is set as 1.0 [units]. In the case of cylindrical domain mean values of the ratio **P(x, t _{i})/P(x, t_{i + 10})** calculated at each point of domain for times i = 0, 10, …, 100 have the average value 1.2699 ± 0.0049. Discrepancy between analytic and numerical solutions computed at the following moments in time equals in average -0.002 ± 0.005 for t = 0, -0.03 ± 0.04 for t = 0.5, -0.05 ± 0.06 for t = 1.0 (see Fig.5c). The analytic solution has been computed on the basis of the equation where the maximum number of used Bessel functions n is 20 (which corresponds to the radius R

Fig. 5. The picture shows a) FEM approximation (linear) of the initial-boundary value problem in the center of the cubic domain i. e. at z = π/2 and at the time t = 0.19 with Δt = 0.01 [units] vs. the analytic result together with b) a distribution of differences between numerical and analytic solutions obtained for each node in Ω domain; c) FEM approximation (linear) of the diffusion equation for the cylindrical domain at z=π/2 and at the following times: t_{0} = 0 (blue), t_{mid} = 0.5 (black) and t_{end} = 1.0 (red) with Δt = 0.01 [units]; d) volume profile of elements within the cylindrical domain where V_{0} = 0.015 and h_{0} = 0.5; mesh quality were enhanced with help of the Metropolis algorithm.

Fig. 6. The picture shows FEM approximation (linear) of the initial-boundary value problem vs. the analytic solution depicted in the center of the cylindrical domain i. e. at z = π/2 ± 0.1 and at the times t = 0 a) and t = 1.0 b) and c). Computations have been done with Δt = 0.01 [units], the diffusion coefficient is set as D = 1.0 [units], and the element size as h_{0} = 0.5 b) and h_{0} = 0.3 c); mesh quality has been enhanced with help of the Metropolis algorithm. To better clarity of the illustration that both solutions (i.e. FEM and anal.) are getting closer each other with decaying h_{0}, in the c) picture the number of presented nodes has been chosen smaller than in the two other cases.