Problem with the fitting

% Define the data
dr_data = [2.5453, 0.042123; 5.0907, 0.075326; 7.636, 0.059506; 10.1813, 0.071553; 12.7267, 0.071365; 15.272, 0.067195; 17.8173, 0.046372; 20.3627, 0.043397; 22.908, 0.017179; 25.4533, -0.0063329; 27.9987, -0.030789; 30.544, -0.047569; 33.0893, -0.089512; 35.6347, -0.080675; 38.18, -0.089138; 40.7253, -0.1102; 43.2707, -0.12061; 45.816, -0.11857; 48.3613, -0.11955; 50.9067, -0.10803; 53.452, -0.10462; 55.9973, -0.099548; 58.5427, -0.097164; 61.088, -0.09994; 63.6333, -0.077017; 66.1787, -0.062839; 68.724, -0.048422; 71.2693, -0.03686; 73.8147, -0.01469; 76.3, 0];
dtheta_data = [2.5453, -0.099251; 5.0907, -0.16064; 7.636, -0.21858; 10.1813, -0.18965; 12.7267, -0.16996; 15.272, -0.18172; 17.8173, -0.15029; 20.3627, -0.12541; 22.908, -0.082786; 25.4533, -0.0071716; 27.9987, 0.03695; 30.544, 0.089002; 33.0893, 0.12873; 35.6347, 0.13092; 38.18, 0.13908; 40.7253, 0.17211; 43.2707, 0.16686; 45.816, 0.15826; 48.3613, 0.14872; 50.9067, 0.15295; 53.452, 0.12677; 55.9973, 0.10964; 58.5427, 0.10223; 61.088, 0.10951; 63.6333, 0.088493; 66.1787, 0.068903; 68.724, 0.054396; 71.2693, 0.035731; 73.8147, 0.030172; 76.3, 0];

% Define initial guess values for parameters
lambda_init = 2.9;
kappa_init = 0.066;
theta_k_init = 0.45;
R_init = 2400;
rout = 76.3;

% Define c
c = sqrt(4 – (rout/R_init)^2);

% Define initial conditions as parameters
a0_bound = 0.0;
b0_bound = 0.0;
c0_bound = 0.0;
d0_bound = 0.0;

% Objective function
objective_function = @(params) compute_error(params, dr_data, dtheta_data, R_init, rout, c);

% Optimization options
options = optimoptions(@fminunc,’Display’,’iter’,’Algorithm’,’quasi-newton’,’MaxFunctionEvaluations’,10000,’MaxIterations’,10000);

% Perform optimization
best_params = fminunc(objective_function, [lambda_init, kappa_init, theta_k_init, a0_bound, b0_bound, c0_bound, d0_bound], options);

% Display best parameters
disp("Optimal Parameters:");
disp("Lambda: " + best_params(1));
disp("Kappa: " + best_params(2));
disp("Theta_k: " + best_params(3));

% Plot the results with the best parameters
[~, dr_fit, dtheta_fit] = compute_error(best_params, dr_data, dtheta_data, R_init, rout, c);

figure;
subplot(2,1,1);
plot(dr_data(:,1), dr_data(:,2), ‘ro’, dr_data(:,1), dr_fit, ‘b-‘);
xlabel(‘r’);
ylabel(‘dr(r)’);
title(‘Fit of dr(r) vs r’);
legend(‘Data’, ‘Fit’);

subplot(2,1,2);
plot(dtheta_data(:,1), dtheta_data(:,2), ‘ro’, dtheta_data(:,1), dtheta_fit, ‘b-‘);
xlabel(‘r’);
ylabel(‘dtheta(r)’);
title(‘Fit of dtheta(r) vs r’);
legend(‘Data’, ‘Fit’);

% Error computation function
function [error, dr_fit, dtheta_fit] = compute_error(params, dr_data, dtheta_data, R, rout, c)
lambda = params(1);
kappa = params(2);
theta_k = params(3);
a0 = params(4);
b0 = params(5);
c0 = params(6);
d0 = params(7);

% Define the function for the differential equations
f = @(r, y) [y(2); ((-1*(lambda+1)*(r*y(2)+y(1)))+(1/r)*((kappa*r^2*cos(theta_k)-(lambda+1))*y(1))-(kappa*r*y(3)*sin(theta_k))+(-16*lambda*r^2)/(c^4*R^2)); y(4); ((-1*(r*y(4)+y(3)))+(1/r)*((kappa*r^2*cos(theta_k)-1)*y(3))+(kappa*r*y(1)*sin(theta_k)))];

% Solve the differential equations using ode45
%r_span = [min(dr_data(:,1)), max(dr_data(:,1))]; % Span of r values from the data
r_span = dr_data(:,1);
[~, sol] = ode45(f, r_span, [a0, b0, c0, d0]);

% Extract solutions
dr_fit = sol(:,1);
dtheta_fit = sol(:,3);

% Compute error (sum of squared differences)
error_dr = sum((dr_data(:,2) – dr_fit).^2);
error_dtheta = sum((dtheta_data(:,2) – dtheta_fit).^2);

% Total error
error = error_dr + error_dtheta;
end

I want to fit the simulated d_r(r) vs r and d_theta(r) vs r values with the above mentioned coupled differential eqns by usuing the fitting parameters lamda, kappa, theta_k. For the sake of good fitting one can use boundary conditions as parameter also. However getting errors. Please help me to solve those errors and fitting those data.
The mathematical eqns and details are below:

$kappa>0$ and $0 leq theta_{k}<frac{pi}{2}$ describe the screening. The boundary conditions are $vec{d}(0)=vec{d}left(r_{text {out }}right)=(0,0)$.

For
$$
c=sqrt{4-left(frac{r_{text {out }}}{R}right)^{2}}
$$

Furthermore,

begin{equation}
(lambda+1)left(r d_{r}^{prime prime}(r)+d_{r}^{prime}(r)right) &+ frac{1}{r}left(kappa r^{2} cos theta_{k}-(lambda+1)right) d_{r}(r)\
& – kappa r d_{theta}(r) sin theta_{k} = -frac{16 lambda r^{2}}{c^{4} R^{2}}
end{equation}

begin{equation}
r d_{theta}^{prime prime}(r) + d_{theta}^{prime}(r) &+ frac{1}{r}left(kappa r^{2} cos theta_{k}-1right) d_{theta}(r)\
& + kappa r d_{r}(r) sin theta_{k} = 0
end{equation}% Define the data
dr_data = [2.5453, 0.042123; 5.0907, 0.075326; 7.636, 0.059506; 10.1813, 0.071553; 12.7267, 0.071365; 15.272, 0.067195; 17.8173, 0.046372; 20.3627, 0.043397; 22.908, 0.017179; 25.4533, -0.0063329; 27.9987, -0.030789; 30.544, -0.047569; 33.0893, -0.089512; 35.6347, -0.080675; 38.18, -0.089138; 40.7253, -0.1102; 43.2707, -0.12061; 45.816, -0.11857; 48.3613, -0.11955; 50.9067, -0.10803; 53.452, -0.10462; 55.9973, -0.099548; 58.5427, -0.097164; 61.088, -0.09994; 63.6333, -0.077017; 66.1787, -0.062839; 68.724, -0.048422; 71.2693, -0.03686; 73.8147, -0.01469; 76.3, 0];
dtheta_data = [2.5453, -0.099251; 5.0907, -0.16064; 7.636, -0.21858; 10.1813, -0.18965; 12.7267, -0.16996; 15.272, -0.18172; 17.8173, -0.15029; 20.3627, -0.12541; 22.908, -0.082786; 25.4533, -0.0071716; 27.9987, 0.03695; 30.544, 0.089002; 33.0893, 0.12873; 35.6347, 0.13092; 38.18, 0.13908; 40.7253, 0.17211; 43.2707, 0.16686; 45.816, 0.15826; 48.3613, 0.14872; 50.9067, 0.15295; 53.452, 0.12677; 55.9973, 0.10964; 58.5427, 0.10223; 61.088, 0.10951; 63.6333, 0.088493; 66.1787, 0.068903; 68.724, 0.054396; 71.2693, 0.035731; 73.8147, 0.030172; 76.3, 0];