Solving odes with ode45 by using actual values for physical constants....

Whenever I try to use actual values for constants in a physical problem I get stuck in a 'busy' message from MATLAB. What do I have to do in order to use the true values for my problem? the code does not work because in inits I`ve got M(r=Rt) = 1e+22 which is a mass of a celestial body. I need to use huge values like that but whenever I do it. I get a bug.

function Europa
% Model for only one layer
%
const=setup;
%mass_constrain = [3.9825 6.372]*1e+21 ;
inits = [0 1e+22];
radious = [1e+3 1];
[r,y] = ode45(@moon,radious,inits,[],const);
rho = const.rho_0*(1+y(:,1)/const.k);
y(:,3) = rho ;
plot(r,y)
legend('P','M','rho')
function europa = moon(r,y,const)
europa = zeros(2,1);
P = y(1) ;
m = y(2) ;
europa(1) = -(const.rho_0*(1+P/const.k)*const.G.*m)./r.^2 ;
europa(2) = 4*pi*(const.rho_0*(1+P/const.k)).*r^2 ;

function const=setup 
const.G = 1e-11 ;
const.rho_0 = 1e+3 ;
const.k = 1e+6 ;