Coupled differential equations in Matlab's ode45

Question

I am trying to numerically solve the following coupled differential equations in Matlab:

$$\begin{array}{l} \frac{df}{dz}=i\delta f(z)+i\kappa b(z)\\ \ \\ \frac{db}{dz}=-i\delta b(z)-i\kappa f(z) \end{array}$$

I have defined these differential equations in a function file:

function [dNdz]= eqs(x,N)
% Parameters:
kappa_D = 0.7;
v = 0.04;
delta = (2*kappa_D/v).*((1./x)-1);
kappa = x.*kappa_D;
% Equations:
dNdz=zeros(2,1);
dNdz(1) = 1i.*delta(x).*N(1)+1i.*kappa(x).*N(2);
dNdz(2) = -1i.*delta(x).*N(2)-1i.*kappa(x).*N(1);
end

Then I tried to solve them using Matlab's ode45 function:

[x,N] =ode45(@(x,N) eqs(x,N),x,[1 0]);

R = abs((N(:,2)./N(:,1))).^2;

plot(x,R)

What I am trying to plot is the absolute-value-squared of the ratio of the two solutions $R$ given above. My results look like this:

But this is not the right solution. I know that from the analytic solutions to these differential equations, R must be:

$$\boxed{R=\frac{\sinh^{2}(\alpha L)}{\cosh^{2}(\alpha L)-\frac{\delta^{2}}{\kappa^{2}}}}$$

This equation produces the plot:

% Parameters
npt = 2^10; x = linspace(0.95,1.05,npt); x = 1./x;
kappa_D = 0.7; v = 0.04;
delta = (2*kappa_D/v).*((1./x)-1);
kappa = x.*kappa_D; a = sqrt(kappa.^2-delta.^2);
L=10;
% Plot
R = ((sinh(a*L)).^2)./(((cosh(a*L)).^2)-((delta.^2)./(kappa.^2)));
plot(x,R);

This is very strange. The numerical solution has to be in agreement with the analytic one. What is the problem here? Is there something wrong with my Matlab syntax or the way I am using ode45?

Any explanation is greatly appreciated.