"Euler method to solve first order differential equation for radionuclei decay"

(*initialisation*)
tau=1.0*10^6;(* second *)
tfinal=10tau;tinit=0;Nstep=1000;
f0=10^6;(*initial value*)

(*define G[t]*)
G[i_]:=-f[i]/tau;
(*end of 
initialisation*)




(*definitions*)
Deltat=(tfinal-tinit)/Nstep;
ilast=Nstep+1;
f[0]=f0;
(*end of definitions*)

t[i_]:=tinit+i*Deltat;
For[i=0,i≤ilast-1,i++,
f[i+1]=f[i]+G[i]Deltat;
];
ttf=Table[{t[i],f[i]},{i,0,ilast-1}];
lptttf=ListPlot[ttf,PlotLabel->"f(t) vs. t",PlotRange->All];

ntheory[t_]:=f0*Exp[-t/tau];
plt=Plot[ntheory[t],{t,tinit,tfinal},PlotRange->All,PlotStyle->Red];
Show[lptttf,plt]