Different results from ode45 and NDSolve?

Question

I am solving nonlinear ODE + integration groups using NDSolve of mma and ode45 of matlab. At last, I found that they give different results, even after I added RungeKutta method, interpolationorder, changed the value of precisiongoal, and accuracygoal of NDSolve. Do the different answer origin from numerical error, or the wrong coding?

The equations I want to solve are as below. I want to solve v(z),qr(z) and qi(z) and plot v(z0)&b(z0)&qi(z0) v.s.z0, where b=Sqrt(1/qr)/k0, z=z0*k0. The initial conditions are, v(0)=1,b0=100,qi(0)=0, N=100.

MMA codes and results (solve v(z),qr(z),qi(z); b=1/qr^2; plot v(z0), b(z0), qi(z0) v.s. z0=z/k0)

Clear["`*"]
np = 8.;
\[Beta] = 9.5*10^-7;
k0 = 2 Pi;
I0I14 = 1.;
b0 = 100.;
z0max = 100000.;
zmax = z0max*k0;
qr0 = 1/(k0*b0)^2;

ek = 6*10^-5*I0I14* v^2*Exp[-2*qr*rho^2];
ep = 1.8*10^-4*1*I0I14^8*(v^2*Exp[-2*qr*rho^2])^8;
er = 1. + ek + ep;
ei = 9.5*10^-7*I0I14^7* (v^2*Exp[-2*qr*rho^2])^7;

{coreA[v_, qr_, qi_, rho_], coreB[v_, qr_, qi_, rho_]} = 
Exp[-qr rho^2] rho RotationMatrix[-qi rho^2].{ei, 1 - er};
int[core_, v_, qr_, qi_?NumericQ] := 
v NIntegrate[core[v, qr, qi, rho], {rho, 0, Infinity}, MaxRecursion -> 40]

intD = \[Beta]/np*I0I14^7;
eqn = With[{v = v@z, qr = qr@z, qi = qi@z}, 
With[{intA = int[coreA, v, qr, qi], intB = int[coreB, v, qr, qi]},
{D[v, z] == -intA qr - intD v^(2*np - 1),
 D[qr, z] == -2 intA qr^2/v - intD qr v^(2*np - 2),
 D[qi, z] == -3 intA qi qr/v + intB qr^2/v - intD qi v^(2*np - 2)
 }]];

bc = {v[0] == 1., qr[0] == qr0, qi[0] == 0.};

sol = NDSolve[{eqn, bc}, {v, qr, qi}, {z, 0., zmax}, InterpolationOrder -> All] 

LogLinearPlot[v[z*k0] /. sol[[2]] // Evaluate, {z, z0min, z0max}]
LogLinearPlot[{Sqrt[1/qr[z*k0]]/k0} /. sol[[2]] // Evaluate, {z, z0min, z0max}]
LogLinearPlot[{qi[z*k0]} /. sol[[2]] // Evaluate, {z, z0min, z0max}]

Matlab codes and results, y(1)=v, y(2)=qr, y(3)=qi:

clear all

k0 = 2 * pi;    nu_0 = 1;                     
b0 = 100;                        
b_0 = b0 * k0;                  
qr_0 = 1/b_0^2;                    % initial qr
qi_0 = 0;                         

options = odeset('RelTol',1e-3,'AbsTol',[1e-4 1e-4 1e-4]);
ts0 = 0; 
tf0 = 1E5;                         
ts = round(0 * k0); 
tf = round(tf0 * k0);
nsteps = 1E5;
step = (tf - ts) / nsteps;
tspan = linspace(ts, tf, nsteps);

T = zeros(nsteps,1);
F = zeros(nsteps,3);
warning('OFF','MATLAB:integral:MinStepSize');
[T,F] = ode45(@test_fun,tspan,y0,options);

%------------Plot--------------%

z0 = T / k0;
nu = F(:,1);
qr = F(:, 2);
qi = F(:, 3);
b0 = 1./(F(:,2)).^0.5 / k0;

plot(z0, nu);      
set(gca,'XScale','log')
plot (z0, b0);     
set(gca,'XScale','log')
plot (z0, qi);    
set(gca,'XScale','log')

-----------function------------
function dy=test_fun(x,y)
I14 = 1;
beta = 9.5 * 10^(-7);
np = 8;

ee = @(r)(I14 * y(1)^2 * exp(-2 * y(2) * r.^2));
ek = @(r) (6 * 10^(-5) * ee(r));
ep = @(r) (1.8 * 10^(-4) * 1 * ee(r).^8);
er = @(r) (1 + ek(r) + ep(r));
ei = @(r) (9.5 * 10^(-7) * ee(r).^7);

% A, B, D
A_func = @(r)(exp(-y(2) * r.^2).* ((1 - er(r)).* sin(y(1) * r.^2 * pi/180) + ei(r).* cos(y(3)* r.^2 * pi/180)).* r);
A = y(1) * integral(A_func, 0, Inf);

B_func = @(r)(exp(-y(2) * r.^2).* ((1 - er(r)).* cos(y(1) * r.^2 * pi/180) - ei(r).* sin(y(3) * r.^2 * pi/180)).* r);
B = y(1) *integral(B_func, 0, Inf);

D = beta * I14^7 / np;


dy = zeros(3,1);
dy(1) = -A * y(2) - D * y(1)^(2 * np - 1);
dy(2) = -2 * A * y(2)^2 / y(1) - D * y(2) * y(1)^(2 * np -2);
dy(3) = -3 * A * y(2) * y(3) / y(1) + B * y(2)^2 / y(1) - D * y(3) * y(1)^(2 * np - 2);

The plots from mma(left) and matlab(right) are: