clc;

% unstable polynom
% P^4+2p^3+3p^2+10p+8

a4=1;
a3=2;
a2=3;
a1=10;
a0=8;

% stable polynom
% P^4+12p^3+24p^2+2p+3
 
% a4=1;
% a3=12;
% a2=24;
% a1=2;
% a0=3;


% Hurwitz determinant
H=[a3 a1 0  0;
   a4 a2 a0 0;
   0  a3 a1 0;
   0  a4 a2 a0];

H1=a3 %det(H1)
H2=[a3 a1;
    a4 a2];
detH2=a3*a2-a4*a1
det(H2)
H3=[a3 a1 0;
   a4 a2 a0;
   0  a3 a1;];
detH3=a3*a2*a1-a3^2*a0-a1^2*a4
det(H3)

detH4=a0*detH3
det(H)








%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Method of derivative order reduction
% a4*y'''' + a3*5y''' + a2*y'' + a1*y' +a0*y = z 
% p3y'=5*p3y +7*p2y -29*py -30*y, y'''(0)=0
% p2y'=p3y, y''(0)=0
% py'=p2y, y'(0)=0
% y'=py, y(0)=0


init=[0;0;0;0];
A=[-a3/a4 -a2/a4 -a1/a4 -a0/a4;
   1 0 0 0;
   0 1 0 0;
   0 0 1 0];

b=[10; 0; 0; 0]; % z... constant signal in the first equation


myodefun = @(t,y) A*y + b;

tspan=[0,10];
tol=1e-5;
ABSTOL(1:length(init))=tol;
options = odeset('RelTol',tol,'AbsTol',ABSTOL);%[1e-10 1e-10 1e-10 1e-10]

[T,Y] = ode23(myodefun,tspan,init, options);

figure
hold on
plot(T,Y(:,4))
hold off
legend('y')
xlabel('t')