# General ODE solved with PHREEQC`s runge-kutta
# The concentration of a chemical element gives the integral.
# The concentration should be about 1 mM for 1e-8 accuracy.

PRINT; -reset false; -user_print true;# -high_precision true
SOLUTION 1
Hdg 1 # use sufficient chemical to balance a negative integral, but not too much for tolerance

RATES
solve_ODE
# integrate y = x*exp(-x) from 0.8 to 5
# x = total_time, dx = time.
# If lower boundary (LB) >= 0, 1st time_step is LB.
30 x = total_time
40 dx = time
50 y = x * exp(-x)
60 if get(2) = 0 then put(1e-3, 2) # adjust factor to add about 1 mmol steps
70 Save y * dx * get(2)
-end

solve_ODE2 # LB < 0
# integrate y = x*exp(-x) from -0.8 to 5
# If LB < 0, 1st time_step = 1e-16 (a very small number), add LB to x in the rate, see below (graph 2)
10 LB = -0.8
20 if get(1) = 0 then put(LB, 1)
30 x = total_time + get(1)
40 dx = time
50 y = x * exp(-x)
60 if get(2) = 0 then put(1e-3, 2) # adjust factor to add about 1 mmol steps
70 save y * dx * get(2)
-end
END

USE solution 1
KINETICS
solve_ODE; -formula Hdg 1; -runge_kutta 6
-time_step 0.8 10*0.042 9*0.42 # for the plot...
# -time_step 0.8 4.2 # gives the same result...
-tol  1e-10
# -bad_step_max 1000
INCREMENTAL_REACTIONS true

USER_GRAPH 1
-axis_titles x y
10 x = total_time
20 if step_no = 1 then put(tot("Hdg"), 3)
30 if step_no = 1 then put((-x - 1) * exp(-x), 4)
40 plot_xy x, x * exp(-x), symbol = None
50 if x < 4.9 then end
60 anal = (-x - 1) * exp(-x) - get(4)
70 print 'Integrate y = x * exp(-x) from x = 0.8 to x = 5'
80 print 'Analytical =', anal, '. PHREEQC Integral =', (tot("Hdg") - get(3)) / get(2), ', difference =', (tot("Hdg") - get(3)) / get(2) - anal
END

USE solution 1
KINETICS 2
solve_ODE2; -formula Hdg 1; -runge_kutta 6
-time_step 1e-20 10*0.08 10*0.5 # for the plot...
# -time_step 1e-20 5.8 # gives the same result...
-tol  1e-10

USER_GRAPH 1; -detach
USER_GRAPH 2
-axis_titles x y
10 x = total_time + get(1)
20 if step_no = 1 then put(tot("Hdg"), 3)
30 if step_no = 1 then put((-x - 1) * exp(-x), 4)
40 plot_xy x, x * exp(-x), symbol = None 
50 if x < 4.9 then end
60 anal = (-x - 1) * exp(-x) - get(4)
70 print EOL$ + 'Integrate y = x * exp(-x) from x = -0.8 to x = 5'
80 print 'Analytical =', anal, '. PHREEQC Integral =', (tot("Hdg") - get(3)) / get(2), ', difference =', (tot("Hdg") - get(3)) / get(2) - anal
END