Hi Mike.
The way you're building phi is creating a matrix, not a vector, since zeros(m,n) creates an mxn matrix.
I think this is what you're trying to do
function [Gamma,phi] = GAMMA(A,B,N)
% A and B are a scalar
phi = A.^(1:N)'; % raises A to powers 1,2,3,...N, creates a column vector
Gamma = diag( B * ones(N,1) ); % create gamma, a diagonal matrix of B
for i=1:N
for j=1:i-1
Gamma(i,j) = B * phi(i-j); % element i-j is A^(i-j)
end
end
end
I tested this with [gam,phi] = GAMMA(2,1,5);
Seems to work alright.