Step response in iterative calculation

I want to plot the step response of this transfer function

Es^2 + As + B / Fs^2 + Cs + D = 90000 / s^2 + 30s + 90000

I use this code,and the output is shown below

H = tf([0 0 90000],[1 30 90000]);
step(H)

In addition,I want to prove this calculation,so I decompose this equation into iterative form

According to this derivation,my code runs thus,and the output is shown below

 t=0.001; %period
 a=0; %y(n)
 b=0; %y(n-1)
 c=0; %y(n-2)

 A=0;
 B=90000;
 C=30;
 D=90000;
 E=0;
 F=1;

 for n = 0:0.01:0.4,
 a = ( heaviside(n)*(E+t*A+t*t*B) + heaviside(n-1)*(-2*E-t*A) + heaviside(n-2)*E + b*(2*F+t*C) -c*F ) / (F+t*C+t*t*D);
 c=b;
 b=a;
 plot(n,a,'*'),
 hold on,
 end

I try to change period,0.01,0.0001,but the curve is still wrong.

How can I correct this problem?