Solving a boundary value problem using bvp4c: forcing a particular form of solution

Hello,

I am looking to replicate some results from a paper, linked here.

The results are a matlab plot that should generate as the figure in the paper, but I am having trouble using bvp5c, as discussed in the paper. The problem is that I cannot get the same form of the equation. Is there a way to force the solution to look more like that in the paper?

I provide some screenshots, and the relevant code I have been using here:

The code I have come up with so far is here:

%% Simulating the Raman Amplifier as done in:
% https://www.osapublishing.org/jlt/abstract.cfm?uri=jlt-33-18-3773
% Based on example of bvp4c

%% Defining the system of differential equations

for i = 2:1:20

    solinit = bvpinit(linspace(0,7e3,i),@(x)init_sol(x));
    options = bvpset('Stats','on','RelTol',1e-9);
    sol = bvp5c(@ode,@bc,solinit,options);
    % The solution at the mesh points
    x = sol.x;
    y = sol.y;

    figure(i);
    clf reset
    plot(x,y','*')
    title('Example problem for MUSN')
    ylabel('bvp4c and MUSN (*) solutions')
    xlabel('x')
    shg

end

function dydx = ode(x,y)
%EX1ODE ODE function for Example 1 of the BVP tutorial.
%   The components of y correspond to the original variables
%   as  y(1) = P_r, y(2) = P_p, y(3) = P_sbs.

g_r = 4.4e-14;
epsilon = 0.55;
g_b = 4e-13;
A_eff = 5.2e-11;
alpha_p = 2e-4;
alpha_r = 3e-4;
lambda_r = 1651e-9;
lambda_p = 1540e-9;

dydx =  [ epsilon*g_r*y(2)*y(1)/A_eff - g_b*y(1)*y(3)/A_eff - alpha_r*y(1) 
          lambda_r/lambda_p * epsilon*g_r*y(2)*y(1)/A_eff + lambda_r/lambda_p * epsilon*g_r*y(2)*y(3)/A_eff + alpha_p*y(2)
          -epsilon*g_r*y(2)*y(3)/A_eff - g_b*y(1)*y(3)/A_eff + alpha_r*y(3) ];

end

function res = bc(ya,yb)
%EX1BC Boundary conditions for Example 1 of the BVP tutorial.
%   RES = EX1BC(YA,YB) returns a column vector RES of the
%   residual in the boundary conditions resulting from the
%   approximations YA and YB to the solution at the ends of 
%   the interval [a b]. The BVP is solved when RES = 0. 
%   The components of y correspond to the original variables
%   as  y(1) = P_r, y(2) = P_p, y(3) = P_sbs.

P_p = 4;
P_r = 6.5e-3;
P_sbs = 1e-6;

res = [ ya(1) - P_r
        yb(2) - P_p
        yb(3) - P_sbs];
    
end

function v = init_sol(x)
%EX1INIT guess for Example 1 of the BVP tutorial.
v = [6.5e-3
             4
             1e-10 ];
end

The desired output should look like this:

Any help or discussion would be greatly appreciated.

[Jaya]

Were you able to solve it later? Because I am also stuck with something similar.