function B = cramermatrix(A,i,b)
%   given matrix A, column nr i and vector b
%   form matrix B as A with column i replaced by b
    B = A;
    B(:,i) = b(:);