I have a function of the form x=x+dx, whereas x(z) and dx(x,z) both depend on a increasing integer. Now I would like to evaluate this function after n time steps. In addition to that I'd like to plot it and compare it to an analytic solution. So, besides the analytic plot there should be points at x0, x1=x0+dx, x2=x1+dx, x3=x2+dx, and so on.

How can this be done? Thank you very much in advance!

A big edit to make it clear what I'd like to achieve.

Given the Saha equation

$\frac{x_e^2}{1-x_e}=\frac{1}{n_b}\left(\frac{2\pi m_ek_BT}{h^2}\right)^{3/2}\exp\left(\frac{-I}{k_BT}\right)$ with $x_e$ the ionisation fraction, $n_B=n_H+n_p=n_e/x_e$ the baryon density (H for hydrogen, e for electrons), $m_e$ the electron mass, $k_B$ the Boltzmann constant, $T$ the temperature, $h$ the Planck constant, and $I=13.6$ eV the binding energy of the ground state of hydrogen.

This equation can be solved by mathematica. Furthermore, we have

$\frac{\mathrm{d}x_e}{\mathrm{d}t}=\alpha\left[-n_Bx_e^2+\left(\frac{2\pi m_ek_BT}{h^2}\right)\left(1-x_e\right)\exp\left(\frac{-I}{k_Bt}\right)\right]$.

On the one hand we have the Saha equation, which can be solved for whatever time, on the other hand we have a numerical change per time step. I wanted to compare a numerical result (solving the Saha equation at a certain time and evolving from there with numerical timesteps) to an analytic result (solving the Saha equation at the different times).

In the end I did it with MATLAB and the following code, where I added a little twist as I calculated it for a certain redshift instead of time steps.

function [ x,y ] = NumSaha()
%% Solve the Saha Equation numerically
% Setting constants in arbitrary units
n0 = 10^-7;
t = 2.7;
j = 13.6*5.291*10^-19;%Planck Units
a = 3*10^-13;
h0 = 1/(4.55*10^17);%inverse of Hubble time

%set initial values
x = zeros(1,2001);
x(2001) = 0.0000626948; %x(z=2000) from mathematica

%numerical solution from x(z=2000) to x(z=0)
while z > 0
nb = n0*(1+z)^3;
x(z) = x(z+1) - a*(-nb*(x(z+1))^2 + (2*pi*m*k*t/h^2)^(3/2)*...
    (1 - x(z+1))*exp(-j/k/t))*(-1/(h0*sqrt(omega)*(1 + z)^2.5));

%now the solution from x(z=0) (given from mathematica / analytical
%solution) to x(z=2000)
y = zeros(1,2001);
y(1) = 0.970117; %x(z=0)

while z < 2001
nb = n0*(1+z)^3;
y(z+1) = y(z) - a*(-nb*(y(z))^2 + (2*pi*m*k*t/h^2)^(3/2)*...
    (1 - y(z))*exp(-j/k/t))*(-1/(h0*sqrt(omega)*(1 + z)^2.5));

%% Analytic Solution
f = @(z) exp(-j/k/t)*(-(2*pi*m*k*t/h^2)^(3/2) +...
    sqrt((2*pi*m*k*t/h^2)^(3) +...

%% Plots
%Plot properties are stored here:

%create the plot
hold on
hold off

%adjust the plot settings
plt = Plot(); % create a Plot object and grab the current figure

plt.XLabel = 'z'; % xlabel
plt.YLabel = 'x_e'; %ylabel
plt.Title = ['Saha Equation']; % plot title

%plt.YScale = 'log';
plt.XScale = 'log';
plt.YLim = [0,0.02];
plt.XLim = [100,2000];

plt.Legend = {'analytical','numerical (increasing)',...
    'numerical (decreasing)'};
plt.LegendBox = 'on';
plt.LegendLoc = 'best';

% %% Save the plot
% plotFile = ['Saha-Equation.png'];
% plt.export(plotFile);

  • $\begingroup$ Can you give us an example analytical form of $x(z)$ and $dx(z)$? As you saw with my last attempt to answer, you had not provided enough information and I misunderstood. I am still unsure exactly of what you need in this case. Are you trying to compare a result obtained recursively to a closed-form analytical expression? Can you provide an actual example of a system you would like to plot? $\endgroup$
    – MarcoB
    Dec 5, 2019 at 16:19
  • 1
    $\begingroup$ To solve a differential equation, you shouldn't be writing your own Euler forward equation like this but just use NDSolve. It will be much more accurate. $\endgroup$ Jan 7, 2021 at 9:46

1 Answer 1


It's a bit unclear what you would like, but try perhaps the following:

    NestList[# + dx &, x, n],
    PlotRange -> {{0, 20}, {0, 310}}
  {{x, 1}, 0, 100},
  {{dx, 0.1}, 0, 10},
  {{n, 2}, 1, 20, 1}
  • $\begingroup$ Thank you very much for your reply! I elaborated a bit more on my problem. :) $\endgroup$
    – kalle
    Dec 4, 2019 at 21:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.