Hello Jaime,
Suppose you want the results up to a polynomial of degree p. The following code creates a (p+1)x(p+1) matrix c of polynomial coefficients, with the coefficent of psi_0 in the first column and the coefficients of psi_p in the (p+1)st column. (Offset by 1 is necessary since Matlab is not zero-based).
The jth column of c has the x^0 coefficient at the bottom, moving upward to the x^j coefficient at the top, with zeros above that. Each column can be used directly in polyval, where the leading zeros could be removed but don't matter.
To check the relation, the columns can be combined directly since the powers of x line up. The coffiicients on the left hand side, representing x*psi(j-1), are shifted up by one to account for the factor of x.
p = 6; % degree of largest polynomial
ain = (1:20)+1; % ain and bin are overlong, doesn't matter
bin = 1:20;
n = 3; % index shift to avoid nonpositive indices
a = [zeros(1,n) ain];
b = [zeros(1,n) bin];
nc = p+n+3; % provide extra space for coefficients matrix
c = zeros(nc,nc);
c(n+1,n) = 1; % psi(0) = 1 including the index shift;
k = n+1;
for j = k:nc
c(k,j) = (1/b(j))*( c(k,j-1) -b(j-1)*c(k-2,j-2) -a(j)*c(k-1,j-1) );
end
for k = n+2:nc
for j = k-1:12
c(k,j) = (1/b(j))*( c(k,j-1) -b(j-1)*c(k-2,j-2) -a(j)*c(k-1,j-1) );
end
end
% trim to put 1 in the upper left corner, then trim to (p+1)*(p+1)
c = c(n+1:end,n:end);
ctemp = c(1:p+1,1:p+1);
c = zeros(p+1,p+1);
% shift the columns downward so that the x^0 term is in the bottom row
for j = 1:ncol
c(p+2-j:end,j) = ctemp(1:j,j);
end
c
% check the relation
% x*psi(j-1) = b(j-1)*psi(j-2) +a(j)*psi(j-1) +b(j)*psi(j)
J = 5; % for example
ofs = 1; % offset since jth column corresponds to psi(j-1)
rhs = bin(J-1)*c(:,J-2+ofs) + ain(J)*c(:,J-1+ofs) +bin(J)*c(:,J+ofs)
lhs = [c(2:end,J-1+ofs);0]
