How do I solve a system using rkf45 method with shooing technique?

I tried to solve a ode system using rk method with shooing technique. But I have a lot of errors. How can I solve this issue.

function bvp_shooting_rk_2()
% Define constants and parameters
rho_f = 783; rCp_f = 2090 * rho_f; k_f = 0.145; rb_f = 21e-6; Pr = 6;
rho_s1 = 1800; rCp_s1 = 717 * rho_s1; k_s1 = 5000; rb_s1 = 6.30e7;
rho_s2 = 10500; rCp_s2 = 235 * rho_s2; k_s2 = 429; rb_s2 = 63e-6;
m = 3; phi_s1 = 0.01; phi_s2 = 0.01; M = 0; K_p = 0.5; beta = 0.5;
S = 0.5; % Initial condition for theta at eta = 0
B_i = 0.5; Q = 0.01; Ec = 0.8;

% Assume values for mu_f and sigma_f, as they are not defined in the problem
mu_f = 1; % Dynamic viscosity of the base fluid (assumed)
sigma_f = 1; % Electrical conductivity of the base fluid (assumed)

rho_hnf = (1-phi_s2)*((1-phi_s1)*rho_f + phi_s1*rho_s1) + phi_s2*rho_s2;
rCp_hnf = (1-phi_s2)*((1-phi_s1)*rCp_f + phi_s1*rCp_s1) + phi_s2*rCp_s2;
mu_hnf = mu_f / ((1-phi_s1)^2.5 * (1-phi_s2)^2.5);
k_bf = k_f * (k_s1 + (m-1)*k_f – (m-1)*phi_s1*(k_f – k_s1)) / (k_s1 + (m-1)*k_f – phi_s1*(k_f – k_s1));
k_hnf = k_bf * (k_s2 + (m-1)*k_bf – (m-1)*phi_s2*(k_bf – k_s2)) / (k_s2 + (m-1)*k_bf – phi_s2*(k_bf – k_s2));
rb_hnf = (1-phi_s2)*((1-phi_s1)*rb_f + phi_s1*rb_s1) + phi_s2*rb_s2;
sigma_bf = sigma_f * (1-phi_s1)^2.5;
sigma_hnf = sigma_bf * (1-phi_s2)^2.5;

A1 = (1-phi_s1)^2.5 * (1-phi_s2)^2.5 * ((1-phi_s2)*(1-phi_s1 + phi_s1*rho_s1/rho_f) + phi_s2*rho_s2/rho_f);
A2 = (1-phi_s2) * (1-phi_s1 + phi_s1*rho_s1/rho_f) + phi_s2*rho_s2/rho_f;
A3 = (1-phi_s2) * (1-phi_s1 + phi_s1*rCp_s1/rCp_f) + phi_s2*rCp_s2/rCp_f;
A4 = (1-phi_s1)^2.5 * (1-phi_s2)^2.5;

% Initial guess for the shooting method
xi_guess = [1; 1]; % Initial guesses for f”(0) and theta'(0)
eta_end = 25; % End of the interval
tol = 1e-6; % Tolerance

% Solve the BVP using the shooting method
[xi, y, eta] = shooting_method(xi_guess, eta_end, tol, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, S);

% Plot the solution
figure;
subplot(2, 1, 1);
plot(eta, y(:, 1), ‘b-‘, ‘LineWidth’, 2);
xlabel(‘eta’);
ylabel(‘f(eta)’);
title(‘Solution for f(eta) using Shooting Method with RK’);
grid on;

subplot(2, 1, 2);
plot(eta, y(:, 4), ‘r-‘, ‘LineWidth’, 2);
xlabel(‘eta’);
ylabel(‘theta(eta)’);
title(‘Solution for theta(eta) using Shooting Method with RK’);
grid on;
end

function [xi, y, eta] = shooting_method(xi_guess, eta_end, tol, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, S)
% Shooting method to adjust xi
options = odeset(‘RelTol’, 1e-6, ‘AbsTol’, 1e-6);

% Initial conditions
y0 = [0; 1; xi_guess(1); S; xi_guess(2)]; % Initial conditions including S for theta at eta = 0

% Solve the ODE system
[eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec), [0 eta_end], y0, options);

% Check the boundary condition at eta_end
bc_error_f_prime = y(end, 2); % f'(eta_end) should be 0
bc_error_theta = y(end, 4); % theta(eta_end) should be 0

% Adjust xi using a simple iteration (could use a more sophisticated method)
xi_new = xi_guess;
iteration_count = 0;
max_iterations = 100; % Limit the number of iterations to prevent infinite loops
while (abs(bc_error_f_prime) > tol || abs(bc_error_theta) > tol) && iteration_count < max_iterations
% Simple adjustment step
xi_new(1) = xi_new(1) – bc_error_f_prime * 0.1;
xi_new(2) = xi_new(2) – bc_error_theta * 0.1;
y0 = [0; 1; xi_new(1); S; xi_new(2)];
[eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, M), [0 eta_end], y0, options);

bc_error_f_prime = y(end, 2);
bc_error_theta = y(end, 4);
iteration_count = iteration_count + 1;
end

if iteration_count >= max_iterations
error(‘Shooting method failed to converge.’);
end

xi = xi_new;
end

function dydeta = ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, M)
% Extract variables from y
f = y(1);
f_prime = y(2);
f_double_prime = y(3);
theta = y(4);
theta_prime = y(5);

% Define the system of ODEs
f_triple_prime = A1 * ((f_prime)^2 – f * f_double_prime) + A2 * M * f_prime + K_p * f_prime – beta * A1 * (f^2 * f_triple_prime – 2 * f * f_prime * f_double_prime);
theta_double_prime = k_f * Pr * (A3 * f * theta_prime – Ec * A4 * (f_double_prime)^2 – Q * theta) / k_hnf;

% Return derivatives as a column vector
dydeta = [f_prime; f_double_prime; f_triple_prime; theta_prime; theta_double_prime];
endI tried to solve a ode system using rk method with shooing technique. But I have a lot of errors. How can I solve this issue.