function testrk(A,b,c)

for p=1:5,
    p
    for l=1:p,  
        for k=0:p-l,
            test=b'*A^k*c.^(l-1);
            fprintf('%g ?= %g\n',test,factorial(l-1)/factorial(l+k))
            if abs(test-factorial(l-1)/factorial(l+k))>eps;
                return
            end
        end
    end
end
