2nd order nonlinear ode

Hello! I'm trying to solve the following problem numerically: .

with the BCs of

y_inf can be any number between 0 and 1. B2, K are paramteres that I'm playing with (it'll be enough to hold the limits of 1e-1 yp to 1e1). mu tilde can be considered to be 1.

I've tried a lot of ways but with no success. I have failed to get convergence in both shooting method and solving a system of equations in a matrix form Ay=b. I tried to solve it with shooting method using implicit Euler. I’ve incorporated the bisection method to find the best initial guess for the derivative (sitting around -0.55). Still, I’m getting some weird stuff. This led me to use a matrix notation and solving system of equations Ay=b. I’m having troubles with convergence here.

I worked here with nodes from 1 to N-1 where the 0th node is y(x=0) and the Nth node is y_inf.

This is my code:

%solve A y = b, where A is a mass matrix
%   y vector of solution
%   b   right-hand side of equation
%   c   upper diagonal of matrix A, c(n) = 0
%   d   diagonal of matrix A
%   a   lower diagonal of matrix A, a(1) = 0

x=linspace(0,50,500);%building a spatial 1D coordinate
deltax=x(end)./(length(x));%delta x
B=[1e-1,5e-1,1e0,1e1,5e1];K=B(end:-1:1);%thermodynamic parameters
alpha=-2.*B./K.*deltax.^2;beta=-1./K.*deltax.^2;
y0=1;%BC at x=0
y_inf=0.7;
max_iter=5;%number of maximum iterations allowed
ysol={zeros(max_iter,length(x))',zeros(max_iter,length(x))',zeros(max_iter,length(x))',zeros(max_iter,length(x))',zeros(max_iter,length(x))'};
%y=zeros(max_iter,length(x))';%defining memory for y
for q=1:length(B)
    B_sec=B(q);K_sec=K(q);alpha_sec=alpha(q);beta_sec=beta(q);
    rho_ref=(1-y_inf).*exp(-sqrt(1/K_sec.*(2.*B_sec+1)).*x)+y_inf;%approximated analytical solution
    ysol{1,q}(:,1)=rho_ref';%starting with an initial guess for y (analytical solution)
    %r=zeros(1,max_iter);%defining memory for r - the spectral radius
    c=[ones(1,length(x)-2),0];%   c   upper diagonal of matrix A, c(n) = 0
    a=[0,ones(1,length(x)-1)];%   a   lower diagonal of matrix A, a(1) = 0
    d=[(-deltax^2-2).*ones(1,length(x)-1),-deltax^2-4];%   d   diagonal of matrix A
    A=diag(d);%creating one's on the main diagonal
    %%%%%creating the mass matrix%%%%%%
     for i=1:length(x)-2
      A(i,i+1)=1;
     end
      for j=1:length(x)-1
       A(j+1,j)=1;
      end
       A(end-1,end)=-1;
       %Linv=inv(tril(A));%inverse lower diagonal matrix
       %H=triu(A);%upper diagonal matrix
       %r=max(abs(eig(Linv*H)));%spectral radius
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%
       b=zeros(max_iter,length(x))';%defining memory for b - right-hand side of equation
       %%%%%creating b%%%%%%
        for k=2:length(x)-2
         b(k,2)=-1*beta_sec*log(rho_ref(k)/y_inf)+(alpha_sec-2)*y_inf+2*rho_ref(end-1);
        end
          b(1,2)=-1*beta_sec*log(rho_ref(1)/y_inf)+(alpha_sec-2)*y_inf+2*rho_ref(end-1)-y0;
          b(end-1,2)=-1*beta_sec*log(rho_ref(end-2)/y_inf)+(alpha_sec-2)*y_inf;
          b(end,2)=-1*beta_sec*log(rho_ref(end-1)/y_inf)+(alpha_sec-3)*y_inf;
          %%%%%%%%%%%%%%%%%%%%%%
          ysol{1,q}(:,2)=A\b(:,2);%solve the problem - 1st iteration
          %r(2)=max(abs(y(:,2)));%spectral radius - 1st iteration

          %%%%%solving the iterative problem%%%%%%
          for l=2:max_iter    
           if infinity_norm(ysol{1,q}(:,l)-ysol{1,q}(:,l-1))>10^(-2) || norm(ysol{1,q}(:,l)-ysol{1,q}(:,l-1))>10^(-2)
            for h=2:length(x)-2
             b(h,l+1)=-1*beta_sec*log(ysol{1,q}(h,l)/y_inf)+(alpha_sec-2)*y_inf+2*ysol{1,q}(end,l);  
            end
              b(1,l+1)=-1*beta_sec*log(ysol{1,q}(1,l)/y_inf)+(alpha_sec-2)*y_inf+2*ysol{1,q}(end,l)-y0;
              b(end-1,l+1)=-1*beta_sec*log(ysol{1,q}(end-1,l)/y_inf)+(alpha_sec-2)*y_inf;
              b(end,l+1)=-1*beta_sec*log(ysol{1,q}(end,l)/y_inf)+(alpha_sec-3)*y_inf;
              ysol{1,q}(:,l+1)=A\b(:,l+1);
              %if y1<y1_inf
              %break
              %end
              %%%%%%%%%%%%%%%%%%%%%%%%%
              %r(l+1)=max(abs(y(:,l+1)));%spectral radius
           else
           end
           %inf_norm(q,l)=infinity_norm(ysol{1,q}(:,l)-ysol{1,q}(:,l-1));
           %two_norm(q,l)=norm(ysol{1,q}(:,l)-ysol{1,q}(:,l-1));
          end
           plot(x,ysol{1,q}(:,end),'linewidth',2);
               ylen = length(ysol{1,q}(:,end));
               Y(1:ylen,q) = ysol{1,q}(:,end);  
               hold on
               lg = ['B = ',num2str(B_sec),'K = ',num2str(K_sec)]; % prepares legend
               lgnd(q,1:length(lg)) = lg;
end
ylabel('$\tilde{\rho}$','Interpreter','latex','Fontsize',20);
xlabel('$\tilde{x}$','Interpreter','latex','Fontsize',20)
legend(lgnd)
axis([x(1) x(end) y_inf y0])


              

I'm trying to solve the problem for different values of B and K but I'm complex numbers (for some reason, y values are becoming negative and this goes to the logarithmic term - not feasible becuase y stands for particle density) also for the simplest case of B=K=1. Iteratively I'm using the fixed point iteration. My spectral radius is 0.81 so I know I should converge but that doesn't happen. rho_ref is a very good initial guess (analytical solution for the lineraized problem). I would be happy to get some help.

Roi