Using this function, I created the graph showing all 4 PDF functions.
function[] = mygraph1(Nmax)
 
% Variables
 
xvals = linspace(0,3,Nmax);
alpha = [0.5,1,1.5,2];
it = 0;
hold on 
 
%iteration
 
while it < 4
    it = it + 1;
    a = alpha(it);
    pdfx = a.*exp(a+xvals-a.*exp(xvals));
    plot(xvals,pdfx)
end
 
hold off
 
end





This gave the following graph:
 
Figure 3 PDF plot of question 2)b)i)
ii)
I made a function that takes a value for alpha and the number of datapoints as an input and then outputs a CDF graph. Using this code and a command line I plotted the 4 CDFs
Here is the function:
function[] = mygraph2(a,Nmax)
 
% Variables
 
xvals = linspace(0,3,Nmax);
it = 0;
cdfxvals = zeros(1,Nmax);
 
% Iteration
 
while it < Nmax
    it = it + 1;
    xval = xvals(it);
    cdfx = integral(@(xval) a.*exp(a+xval-a.*exp(xval)),0,xval);
    cdfxvals(it) = cdfx;
end
 
%Graph Plot
 
plot(xvals,cdfxvals)
 
end

And the command line:
 
To give the graph:
 
Figure 4 Graph showing the CDF for 4 different values of alpha

iii)
For this question I used two functions, mymeans, and mymodes. These functions plot the means or modes respectively against the different values of alpha. 
mymodes:
function[] = mymodes(Nmax)
 
%Variables
 
modevals = zeros(1,Nmax);
avals = linspace(0.1,2,Nmax);
it = 0;
xvals = linspace(0,3,Nmax);
 
%Iterate
 
while it < Nmax
    it = it+1;
    a = avals(it);
    pdf = a.*exp(a+xvals-a.*exp(xvals));
    mode = max(pdf);
    modevals(it)=mode;
end

mymeans:
function[] = mymeans(Nmax)
 
% Variables
 
avals = linspace(0.1,2,Nmax);
it = 0;
meanvals = zeros(1,Nmax);
 
% Iteration

while it < Nmax
    it = it + 1;
    a = avals(it);
    mean = integral(@(x) x.*a.*exp(a+x-a.*exp(x)),0,3);
    meanvals(it)=mean;
end
 
%Graph Plot
 
plot(avals,meanvals)
 
end

Then using a command line, we get a graph showing the means and modes as a function of a.
 
 
Figure 5 A graph showing the mean and mode of the PDF for a varying value of a.
iv) 
Firstly I made the function, mysd, that takes the number of desired iterations as an input and plots a function of standard deviation against a.

function[] = mysd(Nmax)
 
% Variables
 
avals = linspace(0.1,2,Nmax);
it = 0;
sdvals = zeros(1,Nmax);
 
% Iteration

while it < Nmax
    it = it + 1;
    a = avals(it);
    mean = integral(@(x) x.*a.*exp(a+x-a.*exp(x)),0,3);
    sd = sqrt(integral(@(x) ((x-mean).^2).*a.*exp(a+x-a.*exp(x)),0,3));
    sdvals(it) = sd;
end
 
% Plot Graph
 
plot(avals,sdvals)
 
end


And then I made the function myIQR which plots interquartile range as a function of a.


function[]=myIQR(Nmax)

% Variables
 
xvals = linspace(0,3,Nmax);
avals = linspace(0.1,2,Nmax);
iqrvals = zeros(1,Nmax);
it = 0;
 
 
%iteration
 
while it < Nmax
    it = it + 1;
    a = avals(it);
    iqrvals(it) = iqr(a.*exp(a+xvals-a.*exp(xvals)));
end
 
plot(avals,iqrvals)
 
end


And using these two functions I plotted the standard deviation and the interquartile range as a function of a.