Output Estimation Difference Equation

This function is intended to find the output of a system of pulse
transfer function given the input and the system which you can have the
orders (na, nb, d), the previous output and the point at which you want
to estimate the system at. This is like a difference equation
estimation for a system.
The pulse transfer function is written with the shift time operator "z" such as
z*y(t) = y(t+1), z^(-1)*y(t) = y(t-1)

Any discrete system can be written as the given pulse transfer function

z^(-d) B(z) y
G(z) = ------------------------- = -----------
A(z) u
-- A = 1 + a_1 z^-1 + a_2 z^-1 + ... + a_na z^(-na)
-- B = b_0 + b_1 z^-1 + b_2 z^-1 + ... + b_nb z^(-nb)

where

-- y : output of the system
-- u : control action (input to the system)
-- A = 1 + a_1 z^-1 + a_2 z^-1 + ... + a_na z^(-na)
-- B = b_0 + b_1 z^-1 + b_2 z^-1 + ... + b_nb z^(-nb)
-- d : delay in the system

Notice that in difference equation you have "u" inputs as a vector and
"y" may not be estimated yet and hence you can initialize "y" with
zeros and each point estimated should be added to "y" as the present
output depends on the previous outputs and the inputs, then, the
function should be like that

y = zeros(1, length(u))
for m=1:length(u)
y(m) = outputestimation( A, B, d, u, y, m );
end

Function input and output

Inputs:
A = [1, a_1, a_2, a_3, ..., a_na]
B = [b_0, b_1, b _2, b_3, ..., a_nb]
d = number represents the delay
y = column vector that contains the previous outputs (Initially it's zero)
u = column vector that contains the inputs to the system, may be square
wave, sinsoidal function, ...etc.
m = point at which we want to find the output at, it represents time.

Outputs
y : the instantaneous output due to previous ouputs and previous and
present input