# Evaluate and Plot x=x+1

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
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:

%create the plot
figure;
hold on
fplot(f,[0,2000])
plot(x)
plot(y)
hold off

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


• 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? – MarcoB Dec 5 '19 at 16:19

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

Manipulate[
ListPlot[
NestList[# + dx &, x, n],
PlotRange -> {{0, 20}, {0, 310}}
],
{{x, 1}, 0, 100},
{{dx, 0.1}, 0, 10},
{{n, 2}, 1, 20, 1}
]

• Thank you very much for your reply! I elaborated a bit more on my problem. :) – kalle Dec 4 '19 at 21:42