solving system of 4 coupled odes using shooting method and boundary conditions given

Hello! i am trying to solve a boundary value problem with four coupled first order odes, with four initial conditions at r=0 and four boundary conditions at r=10. my code is running but im not getting my desired output and it is not satisfying the conditions. Can someone help me finding where am i going wrong?

function proca_star_shooting_method
    % Clear the workspace and command window
    clear;
    clc;
    format long;

    % Define constants and parameters
    infinity = 10;
    w = 0.817;
    x_init = linspace(1e-5, infinity, 1000);  % More efficient space vector

    % Define initial conditions and boundary conditions
    initial_conditions = [0.394, 0.394, 0, 0];
    boundary_conditions = [1, 0, 0.745, 0];

    % Shooting method with RK4 integration
    options = optimset('TolX', 1e-6, 'Display', 'iter');
    phi_c_shoot = fminbnd(@(phi_c) shooting_error(phi_c, x_init, initial_conditions, boundary_conditions, w), 0, 1, options);

    % Solve the BVP using the optimal phi_c obtained from the shooting method
    [r_data, y_data] = solve_bvp(x_init, initial_conditions, phi_c_shoot, w);

    % Plot the results
    plot_results(r_data, y_data, infinity);
end

function error = shooting_error(phi_c, x_init, initial_conditions, boundary_conditions, w)
    % Solve BVP with given phi_c and calculate the error at the boundary
    [~, y_data] = solve_bvp(x_init, initial_conditions, phi_c, w);
    yb = y_data(:, end);
    error = norm(yb - boundary_conditions');
end

function [r_data, y_data] = solve_bvp(x_init, initial_conditions, phi_c, w)
    % Define the initial conditions for the solver
    y0 = [initial_conditions(1), initial_conditions(2), initial_conditions(3), phi_c];
    
    % Initialize the solution arrays
    r_data = x_init;
    y_data = zeros(4, length(x_init));
    y_data(:, 1) = y0;
    
    h = x_init(2) - x_init(1);  % Step size

    % RK4 integration loop
    for i = 2:length(x_init)
        y_data(:, i) = rk4_step(@(r, y) bsode(r, y, w), r_data(i-1), y_data(:, i-1), h);
    end
end

function y_next = rk4_step(odefun, r, y, h)
    % Perform one step of RK4 integration
    k1 = h * odefun(r, y);
    k2 = h * odefun(r + h/2, y + k1/2);
    k3 = h * odefun(r + h/2, y + k2/2);
    k4 = h * odefun(r + h, y + k3);
    y_next = y + (k1 + 2*k2 + 2*k3 + k4) / 6;
end

function dfdr = bsode(r, y, w)
    % Define the system of ODEs to be solved
    N = 1 - 2 * y(3) / r;
    dfdr = [
        4 * pi * r * 0.745 * (y(4)^2 + y(2)^2 / (N^2 * y(1)^2));
        w * y(4) - 0.745 * y(1)^2 * N * y(4) / w;
        4 * pi * r^2 * (((0.745 * y(1)^2 * N^2 * y(4)^2)^2) / (2 * y(1)^2 * w^2) + 0.5 * 0.745 * (y(4)^2 * N + y(2)^2 / (N * y(1)^2)));
        r^2 * y(2) * w / (y(1)^2 * N^2) - 2 * y(4)
    ];
end

function plot_results(r_data, y_data, infinity)
    % Plot the results of the BVP solution
    figure;
    plot(r_data, y_data(1,:), r_data, y_data(2,:), r_data, y_data(3,:), r_data, y_data(4,:));
    axis([0 infinity -0.5 1.2]);
    title('Functions vs r');
    xlabel('r');
    ylabel('\sigma');
    legend('y1', 'y2', 'y3', 'y4');
    grid on;
end

This is my desired plot. Please help me out!

[J. Alex Lee]

just by counting conditions, if you only have 4 x 1st order ODEs, you can only have 4 BCs, you can't have 8?