function [E_rec,F_rec]= diophantine(A,B,N)
%Solves the Diophantine equation, A and B is modelparameters, N2 is the cost horizon

%Initializing E 
E(1)=1;

%Initializing F
Atilda=[A 0]-[0 A]; %Atilda=A(1-q^-1)
v=zeros(size(Atilda));
v(1)=1;
Aprime=v-Atilda; %Aprime= 1-Atilda
F=Aprime(2:size(Aprime,2));   %F=q(1-Atilda)
m=size(Aprime,2)-size(F,2);
F=[F zeros(1,m)];  %Size F needs to be equal Atilda= size Aprime

%Saving the E & F for each loop, including the initial values
E_rec=zeros(N); %Initialize the E that is recorded
E_rec(1,1)=E;
F_rec(1,:)=F; %Initilize F that is recorded

%Recursion loop for Diophantine
for j=2:N
   r(j)=F(1); 
   R=[E,r(j)]; %updates the new R value, which takes the first element in F and add in
   for i=1:size(F,2)-1 %updates the S polynomial 
      S(i)=F(i+1)-Atilda(i+1)*r(j);
   end
   r=[R zeros(1,N-j)]; % Set R which is E(j+1) so R represent E(j+2) in next loop
   E_rec(j,:)=r;
   E=R;
   F=[S zeros(1,size(Atilda,2)-size(S,2))]; %Set S which is F(j+1) so S represent F(j+2) in next loop

   F_rec(j,:)=F; 
end  

F_rec=F_rec(:,1:size(F_rec,2)-1);   %Erase the last column.In the calculation
%the size of F needs to equal the size of Atilda. This adds a column.

%Check that 1=E*Atilda-q^-1F
%Sum=conv(E,Atilda);
%F=F(1:size(F,2)-1);
%F=[zeros(1,(size(Sum,2)-size(F,2))) F];
%eq=Sum+F