Lowpass Butterworth filter design (own poles' calculations)

I'm trying to calculate coefficients for my Butterworth filter of my own. I've found this article that looks really reliable. But I'm trying to implement it in MATLAB and I'm getting bandstop filter. Please find attached matlab script below:

clc, clear all, close all;
%DIGITAL BUTTERWORTH FILTER DESING  
%Butterworth LP with -3dB at fc and N order and G gain
fs = 10; %frequency sampling
fc = 1;   %cut off frequency
N = 2;      %filter order
G = 1;      %filter gain
wc = 2*pi*fc;


%analog filter design
%poles and zeros   
k = [1:1:N];
%poles searching
%sk = wc * exp(((j*pi)/(2*N))*(2.*k+N-1))
sk = wc * exp(j*(pi/2)) * exp(j*(2.*k+1)*pi/2*N);
z = 0;
G = 1;
[NUM,DEN] = zp2tf(z,sk,G);

%no zeros for LP
Hs = tf(1,DEN)

%filter analysis
figure(1);
subplot(3,1,1)
pzmap(Hs);
subplot(3,1,2)
impulseplot(Hs);
subplot(3,1,3)

stepplot(Hs);    
fvtool(1,DEN);  

I've calculated poles according to equation $10.3.47$ in the article and created transfer function with this poles as a parameters. But fvtool shows bandstop filter instead of lowpass (as described in article). In pzmap I'm getting also pole on right side of $s$-plane so it seems to be unstable. I dont really know how to interpretet this results.

Maks