function [N, exact_N] = radioactive_decay(N_0, tau, dt)

% Program calculates number of nuclei N(t) as a function of time by approximating 
% the solution of differential equation of radioactive decay via the Euler method 

% Declare and explain or variables explicitly:

% Input:
%   N_0    = initial number of nuclei
%   tau    = time constant
%   dt     = time step

time_steps=int32(5*tau/dt); %evolve system in time for up to five decay constants

% Output:
    N(1:time_steps)=0;  % row vector that contains the number of nuclei at each discrete time step
    exact_N(1:time_steps)=0; % N(t) calculated from exact analytical solution of ODE

% Auxiliary variables:
    time(1:time_steps)=0;

N(1)=N_0;  % initial number of nuclei is specified by us when calling function decay
exact_N(1)=N(1);

% Matlab specific shortcut as a replacement for time(i+1)=... in the loop below
%time = linspace(0, 5*tau-dt, 5*tau/dt); 

for i = 1:time_steps-1
    N(i+1) = N(i) - (N(i)/tau) * dt;
    time(i+1)=time(i)+dt;
    exact_N(i) = exact_N(1)*exp(-time(i)/tau);
end


% Plot nuclei vs. time
plot(time, N, '*');
title('Number of Radioactive Nuclei vs. Time');
xlabel('Time(s)');
ylabel('Number of Nuclei');
hold on
plot(time, exact_N);
legend('N(t)','exact N(t)')



% Save results of numerical ODE solution into a file to be processed by
% another program (such as Origin)
data=[time; N];
fid=fopen('decay_numerics.dat','wt');
fprintf(fid,'%6.4f  %12.8f\n',data);
fclose(fid);


% Save values of the exact ODE solution on discrete time grid into a file to 
% be processed by another program (such as Origin)
data=[time; exact_N];
fid=fopen('decay_exact.dat','wt');
fprintf(fid,'%6.4f  %12.8f\n',data);
fclose(fid);

end %radioactive_decay()
