getting singular jacobian error while using bvp4c solver

when i am using a range for phi_c for 0.01 to 0.1 , i am getting a plot. but when im increasing range from 0.01 to 0.8, i am getting error. please help me solve this.

% Define global variable for initial scalar field values
phi_c_values = logspace(log10(0.01), log10(0.1), 100);

% Physical constants
M_pl_squared = 1.0; 
m = 1.0;  % Scalar mass

% Function to define the differential equations
function dydr = bsode(r, y, p)
    G = 1;  % Gravitational constant
    omega = p(1);
    m=1;

    a = y(1);
    alpha = y(2);
    psi0 = y(3);
    phi = y(4);

    dydr = zeros(4, 1);
    dydr(1) = (0.5) * ((a/r) * ((1 - a^2) + 4 * pi * G * r * a * ...
              (psi0^2 * a^2 * (m^2 + omega^2 / alpha^2) + phi^2)));
    dydr(2) = (alpha/2) * (((a^2 - 1)/r) + 4 * pi * r * ...
              (psi0^2 * a^2 * (omega^2 / alpha^2 - m^2) + phi^2));
    dydr(3) = phi;
    dydr(4) = -(1 + a^2 - 4 * pi * r^2 * psi0^2 * a^2 * m^2) * (phi/r) - ...
              (omega^2 / alpha^2 - m^2) * psi0 * a^2;
end

% Function to define boundary conditions
function res = bsbc(ya, yb, p)
    global phi_c;
    res = [ya(1) - 1;         % a(0) = 1
           ya(2) - 1;         % alpha(0) = 1
           ya(3) - phi_c;     % psi0(0) = phi_c
           yb(3);             % psi0(infinity) = 0
           ya(4)];            % phi(0) = 0
end

% Function for calculating ADM mass
function M = ADM_mass(r, a)
    M = (1 - a.^-2) .* (r / 2);
end

% Main loop over phi_c values
omega = 0.864;
infinity = 15;
x_init = linspace(1e-5, infinity, 1000);

% Lists to store maximum masses and radii
max_masses = [];
max_radii = [];
phi_values = [];
mass_values = [];

for phi_c = phi_c_values
    % Initial guess structure
    solinit = bvpinit(x_init, [1, 1, phi_c, 0], omega);
    
    % Solve the BVP
    sol = bvp4c(@bsode, @bsbc, solinit);
    f = sol.y;
    i = find(f(3, :) < 1e-5, 1);
    
    r = sol.x(1:i);
    a = f(1, 1:i);
    
    % Scale alpha and omega
    alpha_last = 1 / f(2, i);
    scaled_alpha = alpha_last * f(2, 1:i);
    scaled_omega = alpha_last * sol.parameters(1);
    
    % Calculate ADM mass
    M = ADM_mass(r, a);
    M_max = 0.633 * M_pl_squared / m;
    M_scaled = M / M_max;
    
    % Find maximum mass and corresponding radius
    max_mass = max(M_scaled);
    target_mass = 0.99 * max_mass;
    index = find(abs(M_scaled - target_mass) == min(abs(M_scaled - target_mass)), 1);
    max_radius = r(index);
    max_M = M_scaled(index);
    
    % Store results
    max_masses(end+1) = max_M;
    max_radii(end+1) = max_radius;
    phi_values(end+1) = phi_c;
    mass_values(end+1) = max_mass;
    
    fprintf('phi_c = %.5f, max_mass = %.5f, max_radius = %.5f\n', phi_c, max_M, max_radius);
end

% Plotting mass-radius curve
figure;
plot(max_radii, max_masses);
title('Scaled ADM Mass (M) vs Radial Coordinate (r)');
xlabel('r (Radius)');
ylabel('M (Scaled Mass)');
grid on;

% Plotting m vs phi
figure;
plot(phi_values, mass_values, 'orange');
title('Scaled ADM Mass (M) vs \phi_c');
xlabel('\phi_c');
ylabel('M (Scaled Mass)');
grid on;

Answer 1

[Torsten]

I can't interprete your solution curves, but it seems to be a problem for the solver that the Scaled Mass goes to 0 with increasing phi_c.

% Define global variable for initial scalar field values
phi_c_values = linspace(0.01,0.15,200);
% Physical constants
M_pl_squared = 1.0; 
m = 1.0;  % Scalar mass
% Main loop over phi_c values
omega = 0.864;
infinity = 15;
x_init = linspace(1e-5, infinity, 1000);
% Lists to store maximum masses and radii
max_masses = [];
max_radii = [];
phi_values = [];
mass_values = [];
solinit = bvpinit(x_init, [1, 1, phi_c_values(1), 0], omega);
options = bvpset('RelTol',1e-8,'AbsTol',1e-8);
for phi_c = phi_c_values    
    
    % Solve the BVP
    sol = bvp4c(@bsode, @(x,y,p)bsbc(x,y,p,phi_c), solinit, options);
    f = sol.y;
    i = find(f(3, :) < 1e-5, 1);
    
    r = sol.x(1:i);
    a = f(1, 1:i);
    
    % Scale alpha and omega
    alpha_last = 1 / f(2, i);
    scaled_alpha = alpha_last * f(2, 1:i);
    scaled_omega = alpha_last * sol.parameters(1);
    
    % Calculate ADM mass
    M = ADM_mass(r, a);
    M_max = 0.633 * M_pl_squared / m;
    M_scaled = M / M_max;
    
    % Find maximum mass and corresponding radius
    max_mass = max(M_scaled);
    target_mass = 0.99 * max_mass;
    index = find(abs(M_scaled - target_mass) == min(abs(M_scaled - target_mass)), 1);
    max_radius = r(index);
    max_M = M_scaled(index);
    
    % Store results
    max_masses(end+1) = max_M;
    max_radii(end+1) = max_radius;
    phi_values(end+1) = phi_c;
    mass_values(end+1) = max_mass;
    
    fprintf('phi_c = %.5f, max_mass = %.5f, max_radius = %.5f\n', phi_c, max_M, max_radius);
    solinit = bvpinit(sol.x,@(x)interp1(sol.x,sol.y.',x),sol.parameters(1));
end

phi_c = 0.01000, max_mass = 0.05268, max_radius = 12.89409
phi_c = 0.01070, max_mass = 0.05977, max_radius = 12.89409
phi_c = 0.01141, max_mass = 0.06718, max_radius = 12.89409
phi_c = 0.01211, max_mass = 0.07490, max_radius = 12.89409
phi_c = 0.01281, max_mass = 0.08291, max_radius = 12.89409
phi_c = 0.01352, max_mass = 0.09116, max_radius = 12.89409
phi_c = 0.01422, max_mass = 0.09965, max_radius = 12.89409
phi_c = 0.01492, max_mass = 0.10817, max_radius = 12.74978
phi_c = 0.01563, max_mass = 0.11703, max_radius = 12.74978
phi_c = 0.01633, max_mass = 0.12605, max_radius = 12.74978
phi_c = 0.01704, max_mass = 0.13520, max_radius = 12.74978
phi_c = 0.01774, max_mass = 0.14446, max_radius = 12.74978
phi_c = 0.01844, max_mass = 0.15381, max_radius = 12.74978
phi_c = 0.01915, max_mass = 0.16323, max_radius = 12.74978
phi_c = 0.01985, max_mass = 0.17269, max_radius = 12.74978
phi_c = 0.02055, max_mass = 0.18218, max_radius = 12.74978
phi_c = 0.02126, max_mass = 0.19169, max_radius = 12.74978
phi_c = 0.02196, max_mass = 0.20118, max_radius = 12.74978
phi_c = 0.02266, max_mass = 0.20967, max_radius = 12.31686
phi_c = 0.02337, max_mass = 0.21909, max_radius = 12.31686
phi_c = 0.02407, max_mass = 0.22846, max_radius = 12.31686
phi_c = 0.02477, max_mass = 0.23777, max_radius = 12.31686
phi_c = 0.02548, max_mass = 0.24700, max_radius = 12.31686
phi_c = 0.02618, max_mass = 0.25615, max_radius = 12.31686
phi_c = 0.02688, max_mass = 0.26520, max_radius = 12.31686
phi_c = 0.02759, max_mass = 0.27415, max_radius = 12.31686
phi_c = 0.02829, max_mass = 0.28298, max_radius = 12.31686
phi_c = 0.02899, max_mass = 0.29169, max_radius = 12.31686
phi_c = 0.02970, max_mass = 0.30028, max_radius = 12.31686
phi_c = 0.03040, max_mass = 0.30874, max_radius = 12.31686
phi_c = 0.03111, max_mass = 0.31706, max_radius = 12.31686
phi_c = 0.03181, max_mass = 0.32523, max_radius = 12.31686
phi_c = 0.03251, max_mass = 0.33189, max_radius = 11.88394
phi_c = 0.03322, max_mass = 0.33979, max_radius = 11.88394
phi_c = 0.03392, max_mass = 0.34755, max_radius = 11.88394
phi_c = 0.03462, max_mass = 0.35516, max_radius = 11.88394
phi_c = 0.03533, max_mass = 0.36262, max_radius = 11.88394
phi_c = 0.03603, max_mass = 0.36992, max_radius = 11.88394
phi_c = 0.03673, max_mass = 0.37708, max_radius = 11.88394
phi_c = 0.03744, max_mass = 0.38331, max_radius = 11.66748
phi_c = 0.03814, max_mass = 0.39018, max_radius = 11.66748
phi_c = 0.03884, max_mass = 0.39690, max_radius = 11.66748
phi_c = 0.03955, max_mass = 0.40347, max_radius = 11.66748
phi_c = 0.04025, max_mass = 0.40989, max_radius = 11.66748
phi_c = 0.04095, max_mass = 0.41537, max_radius = 11.45102
phi_c = 0.04166, max_mass = 0.42152, max_radius = 11.45102
phi_c = 0.04236, max_mass = 0.42753, max_radius = 11.45102
phi_c = 0.04307, max_mass = 0.43340, max_radius = 11.45102
phi_c = 0.04377, max_mass = 0.43913, max_radius = 11.45102
phi_c = 0.04447, max_mass = 0.44473, max_radius = 11.45102
phi_c = 0.04518, max_mass = 0.44898, max_radius = 11.12633
phi_c = 0.04588, max_mass = 0.45434, max_radius = 11.12633
phi_c = 0.04658, max_mass = 0.45958, max_radius = 11.12633
phi_c = 0.04729, max_mass = 0.46470, max_radius = 11.12633
phi_c = 0.04799, max_mass = 0.46969, max_radius = 11.12633
phi_c = 0.04869, max_mass = 0.47456, max_radius = 11.12633
phi_c = 0.04940, max_mass = 0.47807, max_radius = 10.80164
phi_c = 0.05010, max_mass = 0.48274, max_radius = 10.80164
phi_c = 0.05080, max_mass = 0.48731, max_radius = 10.80164
phi_c = 0.05151, max_mass = 0.49176, max_radius = 10.80164
phi_c = 0.05221, max_mass = 0.49610, max_radius = 10.80164
phi_c = 0.05291, max_mass = 0.50034, max_radius = 10.80164
phi_c = 0.05362, max_mass = 0.50323, max_radius = 10.47695
phi_c = 0.05432, max_mass = 0.50731, max_radius = 10.47695
phi_c = 0.05503, max_mass = 0.51129, max_radius = 10.47695
phi_c = 0.05573, max_mass = 0.51452, max_radius = 10.31460
phi_c = 0.05643, max_mass = 0.51833, max_radius = 10.31460
phi_c = 0.05714, max_mass = 0.52205, max_radius = 10.31460
phi_c = 0.05784, max_mass = 0.52503, max_radius = 10.15226
phi_c = 0.05854, max_mass = 0.52860, max_radius = 10.15226
phi_c = 0.05925, max_mass = 0.53208, max_radius = 10.15226
phi_c = 0.05995, max_mass = 0.53448, max_radius = 9.90874
phi_c = 0.06065, max_mass = 0.53783, max_radius = 9.90874
phi_c = 0.06136, max_mass = 0.54110, max_radius = 9.90874
phi_c = 0.06206, max_mass = 0.54429, max_radius = 9.90874
phi_c = 0.06276, max_mass = 0.54640, max_radius = 9.66522
phi_c = 0.06347, max_mass = 0.54947, max_radius = 9.66522
phi_c = 0.06417, max_mass = 0.55247, max_radius = 9.66522
phi_c = 0.06487, max_mass = 0.55540, max_radius = 9.66522
phi_c = 0.06558, max_mass = 0.55725, max_radius = 9.42170
phi_c = 0.06628, max_mass = 0.56007, max_radius = 9.42170
phi_c = 0.06698, max_mass = 0.56282, max_radius = 9.42170
phi_c = 0.06769, max_mass = 0.56551, max_radius = 9.42170
phi_c = 0.06839, max_mass = 0.56712, max_radius = 9.17819
phi_c = 0.06910, max_mass = 0.56971, max_radius = 9.17819
phi_c = 0.06980, max_mass = 0.57224, max_radius = 9.17819
phi_c = 0.07050, max_mass = 0.57471, max_radius = 9.17819
phi_c = 0.07121, max_mass = 0.57609, max_radius = 8.93467
phi_c = 0.07191, max_mass = 0.57847, max_radius = 8.93467
phi_c = 0.07261, max_mass = 0.58079, max_radius = 8.93467
phi_c = 0.07332, max_mass = 0.58254, max_radius = 8.81291
phi_c = 0.07402, max_mass = 0.58476, max_radius = 8.81291
phi_c = 0.07472, max_mass = 0.58613, max_radius = 8.63027
phi_c = 0.07543, max_mass = 0.58827, max_radius = 8.63027
phi_c = 0.07613, max_mass = 0.59036, max_radius = 8.63027
phi_c = 0.07683, max_mass = 0.59157, max_radius = 8.44763
phi_c = 0.07754, max_mass = 0.59358, max_radius = 8.44763
phi_c = 0.07824, max_mass = 0.59552, max_radius = 8.44763
phi_c = 0.07894, max_mass = 0.59660, max_radius = 8.26500
phi_c = 0.07965, max_mass = 0.59847, max_radius = 8.26500
phi_c = 0.08035, max_mass = 0.60029, max_radius = 8.26500
phi_c = 0.08106, max_mass = 0.60121, max_radius = 8.08236
phi_c = 0.08176, max_mass = 0.60296, max_radius = 8.08236
phi_c = 0.08246, max_mass = 0.60465, max_radius = 8.08236
phi_c = 0.08317, max_mass = 0.60544, max_radius = 7.89972
phi_c = 0.08387, max_mass = 0.60706, max_radius = 7.89972
phi_c = 0.08457, max_mass = 0.60864, max_radius = 7.89972
phi_c = 0.08528, max_mass = 0.60928, max_radius = 7.71708
phi_c = 0.08598, max_mass = 0.61079, max_radius = 7.71708
phi_c = 0.08668, max_mass = 0.61224, max_radius = 7.71708
phi_c = 0.08739, max_mass = 0.61365, max_radius = 7.71708
phi_c = 0.08809, max_mass = 0.61414, max_radius = 7.53444
phi_c = 0.08879, max_mass = 0.61548, max_radius = 7.53444
phi_c = 0.08950, max_mass = 0.61655, max_radius = 7.48878
phi_c = 0.09020, max_mass = 0.61711, max_radius = 7.35180
phi_c = 0.09090, max_mass = 0.61834, max_radius = 7.35180
phi_c = 0.09161, max_mass = 0.61953, max_radius = 7.35180
phi_c = 0.09231, max_mass = 0.61997, max_radius = 7.21482
phi_c = 0.09302, max_mass = 0.62108, max_radius = 7.21482
phi_c = 0.09372, max_mass = 0.62214, max_radius = 7.21482
phi_c = 0.09442, max_mass = 0.62247, max_radius = 7.07785
phi_c = 0.09513, max_mass = 0.62346, max_radius = 7.07785
phi_c = 0.09583, max_mass = 0.62369, max_radius = 6.94087
phi_c = 0.09653, max_mass = 0.62461, max_radius = 6.94087
phi_c = 0.09724, max_mass = 0.62549, max_radius = 6.94087
phi_c = 0.09794, max_mass = 0.62559, max_radius = 6.80389
phi_c = 0.09864, max_mass = 0.62640, max_radius = 6.80389
phi_c = 0.09935, max_mass = 0.62716, max_radius = 6.80389
phi_c = 0.10005, max_mass = 0.62714, max_radius = 6.66691
phi_c = 0.10075, max_mass = 0.62783, max_radius = 6.66691
phi_c = 0.10146, max_mass = 0.62848, max_radius = 6.66691
phi_c = 0.10216, max_mass = 0.62833, max_radius = 6.52993
phi_c = 0.10286, max_mass = 0.62890, max_radius = 6.52993
phi_c = 0.10357, max_mass = 0.62943, max_radius = 6.52993
phi_c = 0.10427, max_mass = 0.62916, max_radius = 6.39295
phi_c = 0.10497, max_mass = 0.62961, max_radius = 6.39295
phi_c = 0.10568, max_mass = 0.62963, max_radius = 6.32446
phi_c = 0.10638, max_mass = 0.63000, max_radius = 6.32446
phi_c = 0.10709, max_mass = 0.62973, max_radius = 6.21806
phi_c = 0.10779, max_mass = 0.63002, max_radius = 6.21806
phi_c = 0.10849, max_mass = 0.62964, max_radius = 6.11167
phi_c = 0.10920, max_mass = 0.62986, max_radius = 6.11167
phi_c = 0.10990, max_mass = 0.63002, max_radius = 6.11167
phi_c = 0.11060, max_mass = 0.62951, max_radius = 6.00344
phi_c = 0.11131, max_mass = 0.62960, max_radius = 6.00344
phi_c = 0.11201, max_mass = 0.62915, max_radius = 5.92226
phi_c = 0.11271, max_mass = 0.62915, max_radius = 5.92226
phi_c = 0.11342, max_mass = 0.62862, max_radius = 5.84109
phi_c = 0.11412, max_mass = 0.62852, max_radius = 5.84109
phi_c = 0.11482, max_mass = 0.62783, max_radius = 5.74977
phi_c = 0.11553, max_mass = 0.62765, max_radius = 5.74977
phi_c = 0.11623, max_mass = 0.62742, max_radius = 5.74977
phi_c = 0.11693, max_mass = 0.62659, max_radius = 5.65845
phi_c = 0.11764, max_mass = 0.62626, max_radius = 5.65845
phi_c = 0.11834, max_mass = 0.62533, max_radius = 5.56713
phi_c = 0.11905, max_mass = 0.62490, max_radius = 5.56713
phi_c = 0.11975, max_mass = 0.62386, max_radius = 5.47582
phi_c = 0.12045, max_mass = 0.62334, max_radius = 5.47582
phi_c = 0.12116, max_mass = 0.62276, max_radius = 5.47582
phi_c = 0.12186, max_mass = 0.62156, max_radius = 5.38450
phi_c = 0.12256, max_mass = 0.62087, max_radius = 5.38450
phi_c = 0.12327, max_mass = 0.61956, max_radius = 5.29318
phi_c = 0.12397, max_mass = 0.61876, max_radius = 5.29318
phi_c = 0.12467, max_mass = 0.61790, max_radius = 5.29318
phi_c = 0.12538, max_mass = 0.61641, max_radius = 5.20186
phi_c = 0.12608, max_mass = 0.61543, max_radius = 5.20186
phi_c = 0.12678, max_mass = 0.61381, max_radius = 5.11054
phi_c = 0.12749, max_mass = 0.61271, max_radius = 5.11054
phi_c = 0.12819, max_mass = 0.61124, max_radius = 5.06488
phi_c = 0.12889, max_mass = 0.61000, max_radius = 5.06488
phi_c = 0.12960, max_mass = 0.60825, max_radius = 4.99639
phi_c = 0.13030, max_mass = 0.60686, max_radius = 4.99639
phi_c = 0.13101, max_mass = 0.60496, max_radius = 4.92790
phi_c = 0.13171, max_mass = 0.60341, max_radius = 4.92790
phi_c = 0.13241, max_mass = 0.60135, max_radius = 4.85941
phi_c = 0.13312, max_mass = 0.59963, max_radius = 4.85941
phi_c = 0.13382, max_mass = 0.59739, max_radius = 4.79092
phi_c = 0.13452, max_mass = 0.59549, max_radius = 4.79092
phi_c = 0.13523, max_mass = 0.59306, max_radius = 4.72243
phi_c = 0.13593, max_mass = 0.59096, max_radius = 4.72243
phi_c = 0.13663, max_mass = 0.58874, max_radius = 4.72243
phi_c = 0.13734, max_mass = 0.58598, max_radius = 4.65394
phi_c = 0.13804, max_mass = 0.58352, max_radius = 4.65394
phi_c = 0.13874, max_mass = 0.58051, max_radius = 4.58545
phi_c = 0.13945, max_mass = 0.57777, max_radius = 4.58545
phi_c = 0.14015, max_mass = 0.57446, max_radius = 4.51696
phi_c = 0.14085, max_mass = 0.57141, max_radius = 4.51696
phi_c = 0.14156, max_mass = 0.56818, max_radius = 4.51696
phi_c = 0.14226, max_mass = 0.56433, max_radius = 4.44847
phi_c = 0.14296, max_mass = 0.56068, max_radius = 4.44847
phi_c = 0.14367, max_mass = 0.55637, max_radius = 4.37999
phi_c = 0.14437, max_mass = 0.55220, max_radius = 4.37999
phi_c = 0.14508, max_mass = 0.54770, max_radius = 4.37999
phi_c = 0.14578, max_mass = 0.54250, max_radius = 4.32362
phi_c = 0.14648, max_mass = 0.53718, max_radius = 4.32362
phi_c = 0.14719, max_mass = 0.53131, max_radius = 4.32362
phi_c = 0.14789, max_mass = 0.52443, max_radius = 4.26725
phi_c = 0.14859, max_mass = 0.51691, max_radius = 4.26725
phi_c = 0.14930, max_mass = 0.50826, max_radius = 4.32362
phi_c = 0.15000, max_mass = 0.49673, max_radius = 4.32362

% Plotting mass-radius curve
figure;
plot(max_radii, max_masses);
title('Scaled ADM Mass (M) vs Radial Coordinate (r)');
xlabel('r (Radius)');
ylabel('M (Scaled Mass)');
grid on;

% Plotting m vs phi
figure;
plot(phi_values, mass_values, 'color','red');
title('Scaled ADM Mass (M) vs \phi_c');
xlabel('\phi_c');
ylabel('M (Scaled Mass)');
grid on;

% Function to define the differential equations
function dydr = bsode(r, y, p)
    G = 1;  % Gravitational constant
    omega = p(1);
    m=1;
    a = y(1);
    alpha = y(2);
    psi0 = y(3);
    phi = y(4);
    dydr = zeros(4, 1);
    dydr(1) = (0.5) * ((a/r) * ((1 - a^2) + 4 * pi * G * r * a * ...
              (psi0^2 * a^2 * (m^2 + omega^2 / alpha^2) + phi^2)));
    dydr(2) = (alpha/2) * (((a^2 - 1)/r) + 4 * pi * r * ...
              (psi0^2 * a^2 * (omega^2 / alpha^2 - m^2) + phi^2));
    dydr(3) = phi;
    dydr(4) = -(1 + a^2 - 4 * pi * r^2 * psi0^2 * a^2 * m^2) * (phi/r) - ...
              (omega^2 / alpha^2 - m^2) * psi0 * a^2;
end
% Function to define boundary conditions
function res = bsbc(ya, yb, p, phi_c)
    res = [ya(1) - 1;         % a(0) = 1
           ya(2) - 1;         % alpha(0) = 1
           ya(3) - phi_c;     % psi0(0) = phi_c
           yb(3);             % psi0(infinity) = 0
           ya(4)];            % phi(0) = 0
end
% Function for calculating ADM mass
function M = ADM_mass(r, a)
    M = (1 - a.^-2) .* (r / 2);
end

[T]

is there a different solver that i can use in this case? i want phi_c from 0.01 to 0.8. i want this plot

[Torsten]

What's the corresponding plot in your code above ? How did you get the plot ? Do you have a publication that uses your model ?

The only boundary value solvers in MATLAB are bvp4c and bvp5c.

[T]

Yes I am referring to a paper where this plot is given. Can bvp5c be better than bvp4c?

[Torsten]

Yes I am referring to a paper where this plot is given.

Could you make the plot with the code from above up to the values of phi_c where the code works ? Do the graphs agree up to this point ?

Can bvp5c be better than bvp4c?

Replace the name "bvp4c" by "bvp5c" in your code and you will see. But I doubt that the solver is the reason for your problems.

[T]

Yes it is working till 0.15 value of phi_c and the graph is correct till here.

I tried using bvp5c but not working. I tried changing the tolerance. Perhaps I should use a different way maybe rk4, I don't know.

[Torsten]

Yes it is working till 0.15 value of phi_c and the graph is correct till here.

Can you give the command how to generate the graph up to phi_c = 0.15 ? It's different from the graphs you implemented up to now.

[T]

here!

% Define global variable for initial scalar field values
phi_c_values = logspace(log10(0.01), log10(0.15), 200); % Extended range for phi_c
% Physical constants
M_pl_squared = 1.0; 
m = 1.0;  % Scalar mass

% Main loop over phi_c values
omega = 0.864;
infinity = 15;
x_init = logspace(log10(1e-5), log10(infinity), 2000); % Higher resolution for x_init
% Lists to store maximum masses and radii
max_masses = [];
max_radii = [];
phi_values = [];
mass_values = [];
options = bvpset('RelTol', 1e-10, 'AbsTol', 1e-10); % Higher precision for solver

% Initial solution guess
solinit = bvpinit(x_init, [1, 1, phi_c_values(1), 0], omega);

for phi_c = phi_c_values
    try
        % Solve the BVP
        sol = bvp4c(@bsode, @(x, y, p) bsbc(x, y, p, phi_c), solinit, options);
        f = sol.y;
        
        % Find where the solution sufficiently approaches zero for psi0
        i = find(f(3, :) < 1e-5, 1);
        if isempty(i)
            i = length(f(3, :)); % Use last index if no value falls below threshold
        end
        
        r = sol.x(1:i);
        a = f(1, 1:i);
        
        % Scale alpha and omega
        alpha_last = 1 / f(2, i);
        scaled_alpha = alpha_last * f(2, 1:i);
        scaled_omega = alpha_last * sol.parameters(1);
        
        % Calculate ADM mass
        M = ADM_mass(r, a);
        M_max = 0.633 * M_pl_squared / m;
        M_scaled = M / M_max;
        
        % Find maximum mass and corresponding radius
        max_mass = max(M_scaled);
        target_mass = 0.99 * max_mass;
        index = find(abs(M_scaled - target_mass) == min(abs(M_scaled - target_mass)), 1);
        max_radius = r(index);
        max_M = M_scaled(index);
        
        % Store results
        max_masses(end + 1) = max_M;
        max_radii(end + 1) = max_radius;
        phi_values(end + 1) = phi_c;
        mass_values(end + 1) = max_mass;
        
        fprintf('phi_c = %.5f, max_mass = %.5f, max_radius = %.5f\n', phi_c, max_M, max_radius);
        
        % Update solinit for the next iteration
        solinit = bvpinit(sol.x, @(x) interp1(sol.x, sol.y.', x, 'linear', 'extrap'), sol.parameters(1));
        
    catch ME
        fprintf('Solver failed for phi_c = %.5f: %s\n', phi_c, ME.message);
        break; % Exit loop if solver fails at higher phi_c
    end
end

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 2000 points and the solution are available in the output argument.
The maximum residual is 4.21097e-07, while requested accuracy is 1e-10.

phi_c = 0.01000, max_mass = 0.05267, max_radius = 12.82686
phi_c = 0.01014, max_mass = 0.05399, max_radius = 12.82686
phi_c = 0.01028, max_mass = 0.05537, max_radius = 12.82686
phi_c = 0.01042, max_mass = 0.05679, max_radius = 12.82686
phi_c = 0.01056, max_mass = 0.05824, max_radius = 12.82686
phi_c = 0.01070, max_mass = 0.05972, max_radius = 12.82686
phi_c = 0.01085, max_mass = 0.06124, max_radius = 12.82686
phi_c = 0.01100, max_mass = 0.06280, max_radius = 12.82686
phi_c = 0.01115, max_mass = 0.06439, max_radius = 12.82686
phi_c = 0.01130, max_mass = 0.06601, max_radius = 12.82686
phi_c = 0.01146, max_mass = 0.06768, max_radius = 12.82686
phi_c = 0.01161, max_mass = 0.06938, max_radius = 12.82686
phi_c = 0.01177, max_mass = 0.07112, max_radius = 12.82686
phi_c = 0.01194, max_mass = 0.07290, max_radius = 12.82686
phi_c = 0.01210, max_mass = 0.07471, max_radius = 12.82686
phi_c = 0.01226, max_mass = 0.07657, max_radius = 12.82686
phi_c = 0.01243, max_mass = 0.07847, max_radius = 12.82686
phi_c = 0.01260, max_mass = 0.08041, max_radius = 12.82686
phi_c = 0.01278, max_mass = 0.08240, max_radius = 12.82686
phi_c = 0.01295, max_mass = 0.08442, max_radius = 12.82686
phi_c = 0.01313, max_mass = 0.08650, max_radius = 12.82686
phi_c = 0.01331, max_mass = 0.08861, max_radius = 12.82686
phi_c = 0.01349, max_mass = 0.09077, max_radius = 12.82686
phi_c = 0.01368, max_mass = 0.09297, max_radius = 12.82686
phi_c = 0.01386, max_mass = 0.09522, max_radius = 12.82686
phi_c = 0.01405, max_mass = 0.09752, max_radius = 12.82686
phi_c = 0.01424, max_mass = 0.09976, max_radius = 12.73593
phi_c = 0.01444, max_mass = 0.10215, max_radius = 12.73593
phi_c = 0.01464, max_mass = 0.10459, max_radius = 12.73593
phi_c = 0.01484, max_mass = 0.10708, max_radius = 12.73593
phi_c = 0.01504, max_mass = 0.10961, max_radius = 12.73593
phi_c = 0.01525, max_mass = 0.11220, max_radius = 12.73593
phi_c = 0.01546, max_mass = 0.11484, max_radius = 12.73593
phi_c = 0.01567, max_mass = 0.11752, max_radius = 12.73593
phi_c = 0.01588, max_mass = 0.12026, max_radius = 12.73593
phi_c = 0.01610, max_mass = 0.12305, max_radius = 12.73593
phi_c = 0.01632, max_mass = 0.12589, max_radius = 12.73593
phi_c = 0.01655, max_mass = 0.12879, max_radius = 12.73593
phi_c = 0.01677, max_mass = 0.13174, max_radius = 12.73593
phi_c = 0.01700, max_mass = 0.13474, max_radius = 12.73593
phi_c = 0.01723, max_mass = 0.13779, max_radius = 12.73593
phi_c = 0.01747, max_mass = 0.14090, max_radius = 12.73593
phi_c = 0.01771, max_mass = 0.14391, max_radius = 12.64565
phi_c = 0.01795, max_mass = 0.14712, max_radius = 12.64565
phi_c = 0.01820, max_mass = 0.15039, max_radius = 12.64565
phi_c = 0.01845, max_mass = 0.15371, max_radius = 12.64565
phi_c = 0.01870, max_mass = 0.15708, max_radius = 12.64565
phi_c = 0.01896, max_mass = 0.16051, max_radius = 12.64565
phi_c = 0.01922, max_mass = 0.16400, max_radius = 12.64565
phi_c = 0.01948, max_mass = 0.16753, max_radius = 12.64565
phi_c = 0.01975, max_mass = 0.17112, max_radius = 12.64565
phi_c = 0.02002, max_mass = 0.17477, max_radius = 12.64565
phi_c = 0.02029, max_mass = 0.17847, max_radius = 12.64565
phi_c = 0.02057, max_mass = 0.18204, max_radius = 12.55601
phi_c = 0.02085, max_mass = 0.18584, max_radius = 12.55601
phi_c = 0.02114, max_mass = 0.18970, max_radius = 12.55601
phi_c = 0.02143, max_mass = 0.19361, max_radius = 12.55601
phi_c = 0.02172, max_mass = 0.19757, max_radius = 12.55601
phi_c = 0.02202, max_mass = 0.20158, max_radius = 12.55601
phi_c = 0.02232, max_mass = 0.20564, max_radius = 12.55601
phi_c = 0.02263, max_mass = 0.20975, max_radius = 12.55601
phi_c = 0.02294, max_mass = 0.21370, max_radius = 12.46700
phi_c = 0.02325, max_mass = 0.21791, max_radius = 12.46700
phi_c = 0.02357, max_mass = 0.22217, max_radius = 12.46700
phi_c = 0.02389, max_mass = 0.22647, max_radius = 12.46700
phi_c = 0.02422, max_mass = 0.23082, max_radius = 12.46700
phi_c = 0.02455, max_mass = 0.23521, max_radius = 12.46700
phi_c = 0.02489, max_mass = 0.23965, max_radius = 12.46700
phi_c = 0.02523, max_mass = 0.24390, max_radius = 12.37862
phi_c = 0.02557, max_mass = 0.24842, max_radius = 12.37862
phi_c = 0.02592, max_mass = 0.25299, max_radius = 12.37862
phi_c = 0.02628, max_mass = 0.25759, max_radius = 12.37862
phi_c = 0.02664, max_mass = 0.26222, max_radius = 12.37862
phi_c = 0.02700, max_mass = 0.26665, max_radius = 12.29087
phi_c = 0.02737, max_mass = 0.27136, max_radius = 12.29087
phi_c = 0.02775, max_mass = 0.27611, max_radius = 12.29087
phi_c = 0.02813, max_mass = 0.28088, max_radius = 12.29087
phi_c = 0.02851, max_mass = 0.28568, max_radius = 12.29087
phi_c = 0.02891, max_mass = 0.29026, max_radius = 12.20375
phi_c = 0.02930, max_mass = 0.29512, max_radius = 12.20375
phi_c = 0.02970, max_mass = 0.30000, max_radius = 12.20375
phi_c = 0.03011, max_mass = 0.30491, max_radius = 12.20375
phi_c = 0.03052, max_mass = 0.30984, max_radius = 12.20375
phi_c = 0.03094, max_mass = 0.31452, max_radius = 12.11724
phi_c = 0.03136, max_mass = 0.31949, max_radius = 12.11724
phi_c = 0.03179, max_mass = 0.32448, max_radius = 12.11724
phi_c = 0.03223, max_mass = 0.32919, max_radius = 12.03134
phi_c = 0.03267, max_mass = 0.33421, max_radius = 12.03134
phi_c = 0.03312, max_mass = 0.33923, max_radius = 12.03134
phi_c = 0.03357, max_mass = 0.34426, max_radius = 12.03134
phi_c = 0.03403, max_mass = 0.34901, max_radius = 11.94605
phi_c = 0.03450, max_mass = 0.35405, max_radius = 11.94605
phi_c = 0.03497, max_mass = 0.35909, max_radius = 11.94605
phi_c = 0.03545, max_mass = 0.36384, max_radius = 11.86137
phi_c = 0.03594, max_mass = 0.36888, max_radius = 11.86137
phi_c = 0.03643, max_mass = 0.37392, max_radius = 11.86137
phi_c = 0.03693, max_mass = 0.37866, max_radius = 11.77728
phi_c = 0.03743, max_mass = 0.38369, max_radius = 11.77728
phi_c = 0.03795, max_mass = 0.38870, max_radius = 11.77728
phi_c = 0.03847, max_mass = 0.39341, max_radius = 11.69380
phi_c = 0.03899, max_mass = 0.39840, max_radius = 11.69380
phi_c = 0.03953, max_mass = 0.40308, max_radius = 11.61090
phi_c = 0.04007, max_mass = 0.40805, max_radius = 11.61090
phi_c = 0.04062, max_mass = 0.41268, max_radius = 11.52860
phi_c = 0.04118, max_mass = 0.41761, max_radius = 11.52860
phi_c = 0.04174, max_mass = 0.42221, max_radius = 11.44687
phi_c = 0.04231, max_mass = 0.42709, max_radius = 11.44687
phi_c = 0.04289, max_mass = 0.43195, max_radius = 11.44687
phi_c = 0.04348, max_mass = 0.43647, max_radius = 11.36573
phi_c = 0.04407, max_mass = 0.44097, max_radius = 11.28516
phi_c = 0.04468, max_mass = 0.44574, max_radius = 11.28516
phi_c = 0.04529, max_mass = 0.45018, max_radius = 11.20516
phi_c = 0.04591, max_mass = 0.45489, max_radius = 11.20516
phi_c = 0.04654, max_mass = 0.45927, max_radius = 11.12573
phi_c = 0.04718, max_mass = 0.46392, max_radius = 11.12573
phi_c = 0.04782, max_mass = 0.46822, max_radius = 11.04686
phi_c = 0.04848, max_mass = 0.47249, max_radius = 10.96855
phi_c = 0.04914, max_mass = 0.47704, max_radius = 10.96855
phi_c = 0.04982, max_mass = 0.48123, max_radius = 10.89080
phi_c = 0.05050, max_mass = 0.48540, max_radius = 10.81359
phi_c = 0.05119, max_mass = 0.48952, max_radius = 10.73694
phi_c = 0.05189, max_mass = 0.49390, max_radius = 10.73694
phi_c = 0.05260, max_mass = 0.49795, max_radius = 10.66083
phi_c = 0.05332, max_mass = 0.50195, max_radius = 10.58526
phi_c = 0.05406, max_mass = 0.50592, max_radius = 10.51022
phi_c = 0.05480, max_mass = 0.51014, max_radius = 10.51022
phi_c = 0.05555, max_mass = 0.51402, max_radius = 10.43571
phi_c = 0.05631, max_mass = 0.51786, max_radius = 10.36174
phi_c = 0.05708, max_mass = 0.52165, max_radius = 10.28828
phi_c = 0.05786, max_mass = 0.52541, max_radius = 10.21535
phi_c = 0.05865, max_mass = 0.52912, max_radius = 10.14294
phi_c = 0.05946, max_mass = 0.53279, max_radius = 10.07104
phi_c = 0.06027, max_mass = 0.53641, max_radius = 9.99965
phi_c = 0.06110, max_mass = 0.53999, max_radius = 9.92876
phi_c = 0.06194, max_mass = 0.54352, max_radius = 9.85838
phi_c = 0.06278, max_mass = 0.54701, max_radius = 9.78849
phi_c = 0.06364, max_mass = 0.55046, max_radius = 9.71910
phi_c = 0.06452, max_mass = 0.55357, max_radius = 9.58180
phi_c = 0.06540, max_mass = 0.55693, max_radius = 9.51388
phi_c = 0.06630, max_mass = 0.56024, max_radius = 9.44643
phi_c = 0.06720, max_mass = 0.56350, max_radius = 9.37947
phi_c = 0.06813, max_mass = 0.56671, max_radius = 9.31298
phi_c = 0.06906, max_mass = 0.56959, max_radius = 9.18141
phi_c = 0.07001, max_mass = 0.57271, max_radius = 9.11633
phi_c = 0.07096, max_mass = 0.57578, max_radius = 9.05170
phi_c = 0.07194, max_mass = 0.57851, max_radius = 8.92383
phi_c = 0.07292, max_mass = 0.58148, max_radius = 8.86057
phi_c = 0.07392, max_mass = 0.58411, max_radius = 8.73539
phi_c = 0.07493, max_mass = 0.58697, max_radius = 8.67347
phi_c = 0.07596, max_mass = 0.58957, max_radius = 8.56631
phi_c = 0.07700, max_mass = 0.59218, max_radius = 8.47538
phi_c = 0.07806, max_mass = 0.59474, max_radius = 8.38542
phi_c = 0.07913, max_mass = 0.59723, max_radius = 8.29641
phi_c = 0.08021, max_mass = 0.59967, max_radius = 8.20835
phi_c = 0.08131, max_mass = 0.60204, max_radius = 8.12123
phi_c = 0.08242, max_mass = 0.60434, max_radius = 8.03503
phi_c = 0.08355, max_mass = 0.60658, max_radius = 7.94974
phi_c = 0.08470, max_mass = 0.60875, max_radius = 7.86536
phi_c = 0.08586, max_mass = 0.61065, max_radius = 7.74033
phi_c = 0.08703, max_mass = 0.61267, max_radius = 7.65822
phi_c = 0.08823, max_mass = 0.61461, max_radius = 7.57688
phi_c = 0.08944, max_mass = 0.61628, max_radius = 7.45670
phi_c = 0.09066, max_mass = 0.61806, max_radius = 7.37755
phi_c = 0.09190, max_mass = 0.61975, max_radius = 7.29924
phi_c = 0.09316, max_mass = 0.62115, max_radius = 7.18326
phi_c = 0.09444, max_mass = 0.62265, max_radius = 7.10697
phi_c = 0.09573, max_mass = 0.62385, max_radius = 6.99424
phi_c = 0.09704, max_mass = 0.62515, max_radius = 6.92000
phi_c = 0.09837, max_mass = 0.62613, max_radius = 6.80999
phi_c = 0.09972, max_mass = 0.62720, max_radius = 6.73775
phi_c = 0.10109, max_mass = 0.62795, max_radius = 6.63084
phi_c = 0.10247, max_mass = 0.62876, max_radius = 6.56045
phi_c = 0.10388, max_mass = 0.62924, max_radius = 6.45621
phi_c = 0.10530, max_mass = 0.62958, max_radius = 6.35376
phi_c = 0.10674, max_mass = 0.62996, max_radius = 6.28632
phi_c = 0.10821, max_mass = 0.62999, max_radius = 6.18643
phi_c = 0.10969, max_mass = 0.62984, max_radius = 6.08826
phi_c = 0.11119, max_mass = 0.62952, max_radius = 5.99148
phi_c = 0.11271, max_mass = 0.62900, max_radius = 5.89644
phi_c = 0.11426, max_mass = 0.62839, max_radius = 5.82339
phi_c = 0.11582, max_mass = 0.62749, max_radius = 5.73867
phi_c = 0.11741, max_mass = 0.62631, max_radius = 5.64758
phi_c = 0.11902, max_mass = 0.62488, max_radius = 5.56027
phi_c = 0.12065, max_mass = 0.62316, max_radius = 5.47204

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

phi_c = 0.12230, max_mass = 0.62114, max_radius = 5.38506

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

phi_c = 0.12398, max_mass = 0.61881, max_radius = 5.30202

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

phi_c = 0.12568, max_mass = 0.61597, max_radius = 5.19688

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

phi_c = 0.12740, max_mass = 0.61289, max_radius = 5.11674

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 1.15997e-09, while requested accuracy is 1e-10.

phi_c = 0.12915, max_mass = 0.60935, max_radius = 5.03546

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 1.96208e-09, while requested accuracy is 1e-10.

phi_c = 0.13092, max_mass = 0.60534, max_radius = 4.95769

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 3.41978e-09, while requested accuracy is 1e-10.

phi_c = 0.13271, max_mass = 0.60064, max_radius = 4.85949

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 6.19073e-09, while requested accuracy is 1e-10.

phi_c = 0.13453, max_mass = 0.59541, max_radius = 4.78021

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 1.17598e-08, while requested accuracy is 1e-10.

phi_c = 0.13637, max_mass = 0.58947, max_radius = 4.70422

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 2.37693e-08, while requested accuracy is 1e-10.

phi_c = 0.13824, max_mass = 0.58252, max_radius = 4.60906

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 5.21514e-08, while requested accuracy is 1e-10.

phi_c = 0.14013, max_mass = 0.57468, max_radius = 4.53983

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 1.28017e-07, while requested accuracy is 1e-10.

phi_c = 0.14205, max_mass = 0.56545, max_radius = 4.46186
phi_c = 0.14400, max_mass = 0.55450, max_radius = 4.38907

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

phi_c = 0.14597, max_mass = 0.54116, max_radius = 4.33668

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1810 points and the solution are available in the output argument.
The maximum residual is 1.42595e-08, while requested accuracy is 1e-10.

phi_c = 0.14797, max_mass = 0.52372, max_radius = 4.28868

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1810 points and the solution are available in the output argument.
The maximum residual is 1.03788e-07, while requested accuracy is 1e-10.

phi_c = 0.15000, max_mass = 0.49679, max_radius = 4.33668

% Plotting mass-radius curve
figure;
plot(max_radii, max_masses);
title('Scaled ADM Mass (M) vs Radial Coordinate (r)');
xlabel('r (Radius)');
ylabel('M (Scaled Mass)');
grid on;

% Plotting m vs phi
figure;
plot(phi_values, mass_values, 'color', 'red');
title('Scaled ADM Mass (M) vs \phi_c');
xlabel('\phi_c');
ylabel('M (Scaled Mass)');
grid on;

% Function to define the differential equations
function dydr = bsode(r, y, p)
    G = 1;  % Gravitational constant
    omega = p(1);
    m = 1;
    a = y(1);
    alpha = y(2);
    psi0 = y(3);
    phi = y(4);
    
    dydr = zeros(4, 1);
    dydr(1) = (0.5) * ((a/r) * ((1 - a^2) + 4 * pi * G * r * a * ...
              (psi0^2 * a^2 * (m^2 + omega^2 / alpha^2) + phi^2)));
    dydr(2) = (alpha/2) * (((a^2 - 1)/r) + 4 * pi * r * ...
              (psi0^2 * a^2 * (omega^2 / alpha^2 - m^2) + phi^2));
    dydr(3) = phi;
    dydr(4) = -(1 + a^2 - 4 * pi * r^2 * psi0^2 * a^2 * m^2) * (phi/r) - ...
              (omega^2 / alpha^2 - m^2) * psi0 * a^2;
end

% Function to define boundary conditions
function res = bsbc(ya, yb, p, phi_c)
    res = [ya(1) - 1;         % a(0) = 1
           ya(2) - 1;         % alpha(0) = 1
           ya(3) - phi_c;     % psi0(0) = phi_c
           yb(3);             % psi0(infinity) = 0
           ya(4)];            % phi(0) = 0
end

% Function for calculating ADM mass
function M = ADM_mass(r, a)
    M = (1 - a.^-2) .* (r / 2);
end

[Torsten]

The graphs are different from the one you posted.

[T]

yes a bit different. but this is the best i can get from this solver

[Torsten]

If you get a different solution, you use different equations and/or parameters in my opinion.