You made a mistake on these two lines:
theta = linspace(thetas(k),thetas(k)+thetad,100);
dum = (theta-thetas(k))/thetad;
because although you calculate a range of theta values that depends on the loop iteration, you then subtract thetas(k) from theta, so the values in dum are always going to be the same (they do not change with the loop iteration). The rest of your code then correctly plots all four identical vectors in the same location because that is what you calculated.
Try it yourself, it is easy to print the first few values of dum:
inside the loop, and you will see that you are always calculating exactly the same values (which you then plot, in exactly the same position, so they look like one line).
Anyway, here is a simplified version of your code (edited after Guillaume's comment):
function [Y,X] = BurnFraction_2()
a = 5;
n = 4;
thetas = [-50 -20 0 20];
thetad = 60;
N = numel(thetas);
X = cell(1,N);
Y = cell(1,N);
for k = 1:N
theta = linspace(thetas(k),thetas(k)+thetad,100);
dum = (theta-thetas(k))/thetad;
temp = -a.*dum.^n;
X{k} = theta;
Y{k} = 1-exp(temp);
end
plot(X{1},Y{1},'-',X{2},Y{2},'--',X{3},Y{3},'.',X{4},Y{4},'-.','linewidth',5);
set(gca,'Xlim',[-60 80],'Ylim',[0 1],'fontsize', 18,'linewidth',2);
xlabel('Crank Angle (°)','fontsize', 18);
ylabel('Cumulative Burn Fraction','fontsize', 18);
which makes this:
