function [P,Q] = pade(C,n,m)
% [P,Q] = pade(C,n,m)
% The coefficient lists for Pn(x) and Qm(x) are computed.
% Some combinations of n and m are incompatible.
% C is the list of Maclaurin coefficients, input.
% n is the degree of Pn(x), input.
% m is the degree of Qm(x), input.
% P is the coefficient list for Pn(x), output.
% Q is the coefficient list for Qm(x), output.
C = flipud(C);
A = C(1:n+m+1);
% First solve an m x m system for the Q's
Mq = zeros(m,m);
for j = 1:m,
for k = 1:m,
Mq(j,k) = C(j+k);
end
end
for j = 1:m,
B(j) = - C(n+j+1);
end
if m > 0, Q = Mq\B'; end
Q = [Q;1]';
% Second solve an n+1 x n+1 system for the P's.
for j=1:n+1,
P(j) = A(j);
kmin = max(1,j-m);
for k=j-1:-1:kmin,
P(j) = P(j) + Q(m-(j-k)+1)*A(k);
end
end
P = fliplr(P);