Solving bvp using shooting method to solve coupled first order odes

I was initially trying using python but cant solve using matlab. please help me out. I have attached a sample code and the equations and boundary conditions i need to use to plot mass-radius plot for one star.

scaled equations: dadr1 = (1/2) * a * ((1 - a**2)/r1 + 4 * np.pi * G * r1 *

(psi0_1**2 * a**2 * (1 + 1 / alpha1**2)/(4*np.pi*G) + phi**2))

dalphadr1 = (alpha1/2) * (((a**2 - 1)/r1 + 4 * np.pi * G * r1 *

(psi0_1**2 * a**2 * (1 / alpha1**2 - 1)/(4*np.pi*G) + phi**2)))

dpsidr1 = phi

dphidr1 = -(1 + a**2 - r1**2 * psi0_1**2 * a**2) * (phi/r1) - \

(1 / alpha1**2 - 1)/(4*np.pi*G)**0.5 * psi0_1 * a**2

boundary conditions:ya[0] - 1, # a(0) = 1

ya[1] - 1/omega, # alpha(0) = 1 scaled

ya[2] - phi_c/(4*np.pi*G)**(1/2), # psi0(0) = phi_c scaled

yb[2] - 0, # Allow psi0(inf) to be close to zero, but not exact

yb[3] - 0

onestar2()

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1500 points and the solution are available in the output argument.
The maximum error is 8.25648e-05, while requested accuracy is 1e-06.

The solution was obtained on a mesh of 1500 points.
The maximum error is  8.256e-05. 
There were 77959 calls to the ODE function. 
There were 60 calls to the BC function. 

The parameter is   1.214 
The final omega value is   0.899 

Warning: Ignoring extra legend entries.

Warning: Ignoring extra legend entries.

Warning: Ignoring extra legend entries.

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1500 points and the solution are available in the output argument.
The maximum error is 2.2136e-05, while requested accuracy is 1e-06.

The solution was obtained on a mesh of 1500 points.
The maximum error is  2.214e-05. 
There were 78134 calls to the ODE function. 
There were 60 calls to the BC function. 

The parameter is   1.106 
The final omega value is   0.948 

Warning: Ignoring extra legend entries.

Warning: Ignoring extra legend entries.

Warning: Ignoring extra legend entries.

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1500 points and the solution are available in the output argument.
The maximum error is 1.30051e-05, while requested accuracy is 1e-06.

The solution was obtained on a mesh of 1500 points.
The maximum error is  1.301e-05. 
There were 71957 calls to the ODE function. 
There were 59 calls to the BC function. 

The parameter is   1.095 
The final omega value is   1.066 

Warning: Ignoring extra legend entries.

Warning: Ignoring extra legend entries.

Warning: Ignoring extra legend entries.

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1500 points and the solution are available in the output argument.
The maximum error is 9.98486e-06, while requested accuracy is 1e-06.

The solution was obtained on a mesh of 1500 points.
The maximum error is  9.985e-06. 
There were 76180 calls to the ODE function. 
There were 60 calls to the BC function. 

The parameter is   1.184 
The final omega value is   1.179 

%Program to calculate the inital data
%Considering a spherically symmetric scalar field
%\psi(t,r)=\psi_0(r)exp(i\omega t)
%Line element in polar-areal gauge
%ds^2=-\alpha(r)^2 dt^2+A(r)^2 dr^2+r^2d\Omega^2

function onestar2

clear;
clc;
format long;

infinity = 15;
omega=0.50;
x_init=1e-5:0.01:infinity;
global phi_c
for phi_c=[0.07 0.04 0.02 0.01]

%solinit = bvpinit(linspace(1e-5,infinity,1000),[1 1 phi_c 0],omega);
    solinit = bvpinit(x_init,[1 1 phi_c 0],omega);

    options = bvpset('stats','on','RelTol',1e-6);

    %condition=true;

    %while condition
    
    sol = bvp5c(@bsode,@bsbc,solinit,options);
        %condition = sol.stats.maxerr >= 1e-5;

    %end

    r_data = sol.x;
    f = sol.y;
    
    i=1;
    
    while (f(3,i)>1e-5)
        i=i+1;
    end
    
    r_data = r_data(1:i);

    %Calculating ADM Mass
    A1=f(1,1:i);
    mass=(r_data./2).*(1-A1.^(-2));

    %Scaling alpha and omega
    al_last = 1./f(2,i);
    al_calc = al_last.*f(2,1:i); 
    omega_f = al_last*sol.parameters;

    fprintf('\n');
    fprintf('The parameter is %7.3f \n',...
        sol.parameters)
    fprintf('The final omega value is %7.3f \n',...
        omega_f)
    
    hold on;

    figure(1)
    plot(r_data,mass);
    axis([0 infinity 0 0.7]);
    title('Mass')
    xlabel('r')
    ylabel('M')
    legend('0.07', '0.04', '0.02', '0.01')

    hold on;
    
    figure(2)
    plot(r_data,al_calc);
    axis([0 infinity 0.7 1.1]);
    title('Lapse')
    xlabel('r')
    ylabel('\alpha')
    legend('0.07', '0.04', '0.02', '0.01')

    hold on;
    
    figure(3)
    plot(r_data,f(3,1:i));
    axis([0 infinity 0 0.1]);
    title('ScalarField')
    xlabel('r')
    ylabel('\Psi_0')
    legend('0.07', '0.04', '0.02', '0.01')

hold off


% figure(4)
% plot(r_data,f(4,:));
% axis([0 infinity 0 2]);
% title('FieldDerivative')
% xlabel('r')
% ylabel('Q')

end
end
% --------------------------------------------------------------------------
%Function to decalre the complete set of differential equations
%y=[y(1),y(2),y(3),y(4)]=[A,alpha,psi_0,Q]
%V=1/2 m^2 \psi\psi* (free field term)
%dV/d\psi=1/2 m^2 \psi*


function dfdr = bsode(r,y,p)
   
m = 1;
%omega = 0.97;
dfdr = [ (1/2)*((y(1)/r)*(1-y(1)^2)+4*pi*r*y(1)*(y(3)^2*y(1)^2*(m^2+p(1)^2/y(2)^2)+y(4)^2))
           (y(2)/2)*((y(1)^2-1)/r+4*pi*r*(y(3)^2*y(1)^2*(p(1)^2/y(2)^2-m^2)+y(4)^2))
           y(4)
           -(1+y(1)^2-4*pi*r^2*y(3)^2*y(1)^2*m^2)*(y(4)/r)-(p(1)^2/y(2)^2-m^2)*y(3)*y(1)^2];

% --------------------------------------------------------------------------
%Function file for declaring the boundary conditions
%Regularity at origin : A(0)=1, \psi_0(0)=\psi_c, Q(0)=0
%Asymptotic flat condition: A(r->\inf)=1, \alpha(r->inf)=1, \psi_0(r->inf)=0 
end
function res = bsbc(ya,yb,p)

global phi_c

res = [ya(1)-1 %fixed
       %yb(1)-1 
       ya(2)-1 %fixed
       %yb(2)-1 
       ya(3)-phi_c %fixed
       yb(3)    %fixed
       ya(4) ]; %fixed


end

[Torsten]

Technically, your code works.

We cannot tell what makes it hard for "bvp5c" to solve your ODE system since we don't know where your equations stem from.

[T]

I have scaled the equations and wrote this code keeping in mind the above code. Can you tell what am i doing wrong?

%Program to calculate the inital data
%Considering a spherically symmetric scalar field
%\psi(t,r)=\psi_0(r)exp(i\omega t)
%Line element in polar-areal gauge
%ds^2=-\alpha(r)^2 dt^2+A(r)^2 dr^2+r^2d\Omega^2
onestar2()

Error using bvp4c (line 196)
Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in solution>onestar2 (line 30)
sol = bvp4c(@bsode,@bsbc,solinit,options);

function onestar2

clear;
clc;
format long;

infinity = 15;
omega=0.50;
x_init=1e-5:0.01:infinity;
global phi_c
for phi_c=[0.07 0.04 0.02 0.01]

%solinit = bvpinit(linspace(1e-5,infinity,1000),[1 1 phi_c 0],omega);
    solinit = bvpinit(x_init,[1 1 phi_c 0],omega);

    options = bvpset('stats','on','RelTol',1e-6);

    %condition=true;

    %while condition
    
    sol = bvp4c(@bsode,@bsbc,solinit,options);
        %condition = sol.stats.maxerr >= 1e-5;

    %end

    r_data = sol.x;
    f = sol.y;
    
    i=1;
    
    while (f(3,i)>1e-5)
        i=i+1;
    end
    
    r_data = r_data(1:i);

    %Calculating ADM Mass
    A1=f(1,1:i);
    mass=(r_data./2).*(1-A1.^(-2));

    %Scaling alpha and omega
    al_last = 1./f(2,i);
    al_calc = al_last.*f(2,1:i); 
    omega_f = al_last*sol.parameters;

    fprintf('\n');
    fprintf('The parameter is %7.3f \n',...
        sol.parameters)
    fprintf('The final omega value is %7.3f \n',...
        omega_f)
    
    hold on;

    figure(1)
    plot(r_data,mass);
    axis([0 infinity 0 0.7]);
    title('Mass')
    xlabel('r')
    ylabel('M')
    legend('0.07', '0.04', '0.02', '0.01')

    hold on;
    
    figure(2)
    plot(r_data,al_calc);
    axis([0 infinity 0.7 1.1]);
    title('Lapse')
    xlabel('r')
    ylabel('\alpha')
    legend('0.07', '0.04', '0.02', '0.01')

    hold on;
    
    figure(3)
    plot(r_data,f(3,1:i));
    axis([0 infinity 0 0.1]);
    title('ScalarField')
    xlabel('r')
    ylabel('\Psi_0')
    legend('0.07', '0.04', '0.02', '0.01')
hold off



% figure(4)
% plot(r_data,f(4,:));
% axis([0 infinity 0 2]);
% title('FieldDerivative')
% xlabel('r')
% ylabel('Q')

end
end
% --------------------------------------------------------------------------
%Function to decalre the complete set of differential equations
%y=[y(1),y(2),y(3),y(4)]=[A,alpha,psi_0,Q]
%V=1/2 m^2 \psi\psi* (free field term)
%dV/d\psi=1/2 m^2 \psi*


function dfdr = bsode(r1,y,p)
   
m = 1;
G = 1;
%omega = 0.97;
a=y(1);
alpha1 = y(2);
psi0_1 = y(3);
phi = y(4);
dadr1 = (1/2) * a * ((1 - a^2)/r1 + 4 * pi * G * r1 * (psi0_1^2 * a^2 * (1 + 1 / alpha1^2)/(4*pi*G) + phi^2));
dalphadr1 = (alpha1/2) * (((a^2 - 1)/r1 + 4 * pi * G * r1 * (psi0_1^2 * a^2 * (1 / alpha1^2 - 1)/(4*pi*G) + phi^2)));
dpsidr1 = phi;
dphidr1 = -(1 + a^2 - r1^2 * psi0_1^2 * a^2) * (phi/r1) - (1 / alpha1^2 - 1)/(4*pi*G)^0.5 * psi0_1 * a^2;
dfdr = [dadr1;dalphadr1;dpsidr1;dphidr1];
%a , alpha1, psi0_1, phi = y
%dadr1 = (1/2) * a * ((1 - a^2)/r1 + 4 * np.pi * G * r1 * (psi0_1^2 * a^2 * (1 + 1 / alpha1^2)/(4*np.pi*G) + phi^2))
%dalphadr1 = (alpha1/2) * (((a^2 - 1)/r1 + 4 * np.pi * G * r1 * (psi0_1^2 * a^2 * (1 / alpha1^2 - 1)/(4*np.pi*G) + phi^2)))
%dpsidr1 = phi
%dphidr1 = -(1 + a^2 - r1^2 * psi0_1^2 * a^2) * (phi/r1) - (1 / alpha1^2 - 1)/(4*np.pi*G)^0.5 * psi0_1 * a^2
end


% --------------------------------------------------------------------------
%Function file for declaring the boundary conditions
%Regularity at origin : A(0)=1, \psi_0(0)=\psi_c, Q(0)=0
%Asymptotic flat condition: A(r->\inf)=1, \alpha(r->inf)=1, \psi_0(r->inf)=0 

function res = bsbc(ya,yb,p)

global phi_c

res = [ya(1)-1 %fixed
       %yb(1)-1 
       ya(2)-1 %fixed
       %yb(2)-1 
       ya(3)-phi_c %fixed
       yb(3)    %fixed
       ya(4) ]; %fixed

end

[Torsten]

I changed some settings in your code:

np.pi is pi in MATLAB,

a , alpha1, psi0_1, phi = y is not a valid assignment

function outputs must be defined: dfdr = [dadr1,dalphadr1;dpsidr1;dphidr1];

functions must end with an "end".

I ended up with an undefined variable G.

[T]

G can be taken as 1. i am getting other errors as well.

[Torsten]

I removed the syntax errors, but MATLAB is not able to solve your equations (see above).

Since you fixed the values of all variables at x = 0, you should try an initial-value solver like ode45 to see what you get for different values of p. But I doubt that fixing all variables at x = 0 is the correct choice of the boundary conditions because it seems you have a singular point there.

[T]

the code i provided in the question is for mass-radius plot for a single boson star for different phi_c values. i want to plot the mass-radius plot for a system of boson stars for different range of phi_c values where each point of the plot denotes one star with one phi_c value. can i do that? if yes, how. i want a single curve which represents the mass-radius plot for system of stars. the final should resemble the blue line. the code should work in such a way that the odes are being solved with respect to the boundary conditions over the range of the phi_c given and finally it plots the final maximum mass and corresponding radius obtained for each phi_c value after which the graph gets stable for one star and plot it in one curve.