Define a step input, and fit OE model as this model structure fits your data the best. See updated code below.
% generate data: sum of 2 exp
d = 10; % input delay
p = 100; % start step
N = 1000;
x = [zeros(p+d,1); (1:(N-p-d))'];
y0 = exp(-0.01*x)+0.8*exp(-0.1*x); % positive sum of two decaying exp
y = y0+0.01*max(y0)*randn(N,1); % add a bit of noise
y = y - y0(1); % remove start value, as fit model around steady state zero
u = [zeros(p,1); ones(N-p,1)];
yid = iddata(y,u,1);
% apply SYSID tools: oe
n = 2;
[sys1,ic1] = oe(yid,[n n d]);
[yh1, FIT1] = compare(yid,sys1);
% apply SYSID tools: oe with initial estimate
init_sys = idpoly(zpk(0.8,[0.99 0.9],-y0(1),1)); % initial sys
init_sys.nk = d; % number of samples input delay
opt = oeOptions('InitialCondition','zero'); % initial state
[sys2,ic2] = oe(yid,init_sys,opt);
[yh2, FIT2] = compare(yid,sys2);
% show results
figure(1);
clf;
t = 0:N-1;
plot(t,u,t,y+y0(1),t,yh1.y+y0(1),t,yh2.y+y0(1),'LineWidth',4);
legend('input','data','SYSID1','SYSID2');
disp('system poles');
disp(roots(sys2.f));
