Diophantine Equation Solver

This function is used to solve the Diophantine euqation

This function is intended to solve the Diphantine equation in the form of
AR + z^(-d) BS = A0Am = alpha;
where

-- 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)
-- R = 1 + r_1 z^-1 + r_2 z^-1 + ... + r_nr z^(-nr)
-- S = s_0 + s_1 z^-1 + s_2 z^-1 + ... + s_ns z^(-ns)
-- d : delay in the system. Notice that this form of the Diaphontaing solution
is available for systems with d>=1

-- alpha = 1 + alpha1 z^-1 + alpha2 z^-1 + ... + alpha_(nalpha z)^(-nalpha) = Am*A0, required characteristic polynomial
-- Am = required polynomial of the model;
-- A0 = observer polynomail for compensation of the order

The function input outputs are given in the following

function [ S, R ] = Diophantine( A, B, d, alpha )

Inputs
A = [1, a_1, a_2, a_3, ..., a_na]
B = [b_0, b_1, b _2, b_3, ..., a_nb]
d = delay time, a number.
alpha = [1, alpha_1, alpha _2, alpha_3, ..., alpha_nalpha], nalpha is
the final order of the closed loop transfer function

Outputs
S = [s_0, s_1, s _2, s_3, ..., s_ns]
R = [1, r_1, r_2, r_3, ..., r_nr]

to find the oreders of the polynomials we use these equations

nr = nb + d - 1
ns = na - 1
nalpha = na + nb + d - 1

the functions is used to estimate the polynomials S and R which are the
numerator and the denomenator of the controller transfer function,
respectively.

The Solution is given in matrix form by solving a linear system of
equations such as

M*theta = (V-Y) --> theta = M^(-1)*(V-Y)

-- M : Sylvester matrix
-- V: vector contains the "alpha" polynomail coefficients without "1" at the
first of it.

V = transpose([alpha_1, alpha _2, alpha_3, ..., alpha_nalpha])
size(nalpha, 1)

-- Y: vector contains the "A" polynomail coefficients without "1" at the
first of it.

Y = transpose([a_1, a_2, a_3, ..., a_na, 0, 0, ..., 0])
size(nalpha, 1)

-- theta: vector contains the unknowns. That is, the coefficients of the R
polynomial and the coefficients of the S polynomial

theta = tranpose([r_1, r_2, ..., r_nr, s_0, s1, s2, ..., s_ns])

An example is added to illustrate how to use the function