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

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

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
k=1;
h=2*pi;
m=1/2000;
a = 3*10^-13;
h0 = 1/(4.55*10^17);%inverse of Hubble time
omega=0.3;

%set initial values
z=2000;
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));
z=z-1;
end

%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)
z=1;

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));
z=z+1;
end

%% 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) +...
    4*exp(j/k/t)*n0*(1+z)^3*(2*pi*m*k*t/h^2)^(3/2)))/(2*n0*(1+z)^3);

%% Plots
%Plot properties are stored here:
addpath('.\lib')

%create the plot
figure;
hold on
fplot(f,[0,2000])
plot(x)
plot(y)
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];

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

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

end

Thank you for opening this question. I think it is a very simple one, unfortunately I am not that experienced using mathematica, nevertheless a simple solution is highly appreciated. :)

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!

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!