Plotting a simple orbit

Hi, I am having a little trouble setting up a simple orbit code so I could plot the orbit. Could anyone see my code? I have all the constants inside "K" and mi is the reduced mass of two bodies, gamma is Gm1m2. The problem is that I should get "r" constant since it's a circular orbit. My orbit equation is written as a comment down there and u = 1/r so the equation is simpler. I set up broke up the 2nd order diferential equation into y1 and y2 and so on... I do not see any mistakes though my r is not constant as well as the velocity (Jupiter's).

function SunsOrbit
% Simple 2 body problem with Jupiter and the Sun 
% CIRCULAR ORBIT AND JUPITER's VELOCITY CONSTANT WERE SUPPOSED
%
%Setting up Constants
K.d=7.785e8 ; %km distance Jupiter-Sun
K.Mj=1.898e27; %kg
K.Msun= 1.989e30; %kg
K.mi=K.Mj*K.Msun/(K.Mj+K.Msun);
K.Vj=47052; %km/h Jupiter's velocity
K.Pj=K.Mj*K.Vj; % Jupiter's linear momentum
K.l=K.d*K.Pj; % Angular Momentum Modulus
K.G=6.67384e-11; % m3 kg-1 s-2
K.gamma = -K.G*K.Mj*K.Msun;
%
Inits=[K.d,K.Vj];
ThetaRange=linspace(0,pi,1000);
[Theta,y]=ode45(@orbit,ThetaRange,Inits,[],K); % y(1) is 1/r(theta)
r=-2*y(:,2)./y(:,1); % THIS SHOULD BE CONSTANT!! 
Vj=-y(:,1).*r;
%[u,i]=pol2cart(Theta,Vj);
%plot(r,Vj);

% Orbit's equation

function da=orbit(Theta,y,K)
da=[y(2);-y(1)+K.mi*K.gamma/K.l^2]; % u'' = -u + mi*gamma/l 
                                    %