Pade Approximation Script - function P not being output correctly

function [a,b,q,p] = padeapp(f,m,n) N = m+n; syms x

for i = 1:N+1 temp = (diff(f,i-1)); x= 0; a(i) = eval(temp)/factorial(i-1); end p(1) = a(1); q(1) = 1; for i = 1:N for j = 1:i if j<= n b(i,j) =0; end if i<=n b(i,i) = 1; end end for j = i+1:N b(i,j) = 0; end for j = 1:i if j<=m b(i,n+j) = -a(i+1-j); end end for j = n+i+1:N b(i,j) = 0; end

    b(i,N+1) = a(i);
end
for i= n+1:N-1
    c = b(:,i);

    c2 = max(c);
    k = find(c == c2);

    if b(k,i) == 0
        fprintf('system is singular,1 i equals %d and k equals %d \n',i,k);
    end
    if k ~= i
        for j = i+1:N+1
            Btemp = b(i,j);
            b(i,j) = b(k,j);
            b(k,j) = Btemp;
        end
    end

    for j = i+1:N
        xm = b(j,i)/b(i,i);
        for k = i+1:N+1
            b(j,k) = b(j,k) -xm*b(i,k);
            b(j,i) = 0;
        end
    end

    if b(N,N) == 0
        fprintf('System is singular,2')
    end

    if m>0
        q(m) = b(N,N+1)/b(N,N);
    end

    for i = N-1:n+1
        tempbq = 0;

        for j = i+1:N
            tempbq = tempbq + b(i,j)*q(j-n);
        end
        q(i-n) = (b(i,N+1) - tempbq)/b(i,i);

        for j = 1:n
            p(n+1-j)= b(n+1-j,N+1)-tempbq;
        end
    end

end

OUTPUT:

[a,b,q,p] = padeapp(sin(x),3,3)

a =

Columns 1 through 6

         0    1.0000         0   -0.1667         0    0.0083

Column 7

         0

b =

Columns 1 through 6

    1.0000         0         0         0         0         0
         0    1.0000         0   -1.0000         0         0
         0         0    1.0000         0   -1.0000         0
         0         0         0    0.1667         0   -1.0000
         0         0         0         0    0.1667         0
         0         0         0         0         0    0.1167

Column 7

         0
    1.0000
         0
   -0.1667
         0
         0

q =

     1     0     0

p =

     0

[Victor]

What am I doing wrong to get this P output? is Q right?