coupled ode for 2nd order

I have coupled differential equations:

a(r) \sin ^{4} f+2 \cot (f) a_{r} f_{r}+\frac{a_{r}}{r}-a_{r r}=0,{l}\lambda_{0} r(1-\cos f) \sin f\left(\cos f-\cos ^{2} f+\sin ^{2} f\right) \\ \quad+\frac{a^{2} \cos f \sin f}{r}-\frac{\cot f \csc ^{2} f a_{r}^{2}}{r}-f_{r}-r f_{r r}=0

with certain boundary conditions.

[Torsten]

Please put your equations, initial and boundary conditions into a readable format.

Answer 1

[Torsten]

I didn't compare equations and boundary conditions with those listed above.

% Defining parameters

delta = 0.02;        % Lower integral bound
R = 5;             % Upper integral bound
theta = 0;          % ArcTan(q/g)
maxPoints = 1e6;    % Maximum numer of grid point used by bvpc4
initialPoints = 10; % Number of initial grid points used by bvpc4
tol = 1e-3;         % Maximum allowed relative error.
L = 10;
N = 2;
n = 0;
m = 0;
g = 5;
lambda = 0;


% Boundary conditions
y0 = [0, -1, N*pi, 0];


% Initial conditions 

A = @(xi) (1-tanh(((L*xi)/R)-(L/3)))/2;
dA = cosh(theta)*(coth(delta)-delta*csch(delta)^2);
F = @(xi) (1+tanh(((L*xi)/R)-(L/3)))/2;
dF = (1-delta*coth(delta))*csch(delta);
solinit = bvpinit(linspace(delta, R, initialPoints), [A(delta), F(delta), dA,dF]);


% Solves system using bvpc4
options = bvpset('RelTol', tol, 'NMax', maxPoints); % This function sets the allowed
    %relative error and maximum number of grid points.  
sol = bvp4c(@(xi, y) heatGauge(xi, y, lambda, g, m, n), @(ya, yb) bcheatGauge(ya, yb, y0),...
    solinit, options);
xi = linspace(delta, R, 1e4);
y = deval(sol, xi);
plot(xi,y)

function dy1 = heatGauge(xi, y, lambda, g, m, n)
  dy1 = [y(3)...
         y(4)...
         y(3)./xi + (g^2 * (1+y(1)) * (1+(lambda^2*y(4)^2)) * sin(y(2))^2)...
         (1./(1+(lambda^2*(n*(y(1)+1)./xi).^2*sin(y(2))^2))) .* ( ((sin(y(2))*cos(y(2))*(n*(y(1)+1)./xi).^2) + (m^2*sin(y(2)))) - (y(4)./xi).*( ((lambda^2*(n*(y(1)+1)./xi)*sin(y(2))^2).*((n*(y(1)+1)./xi)+(2*xi.*(((xi*n*y(3))-(n*y(1))-n)./xi.^2)))) + 1 + (lambda^2*y(4)*xi.*(n*(y(1)+1)./xi).^2.*sin(y(2))*cos(y(2))) ) ) ];
end

function res = bcheatGauge(ya, yb, y0)
  res = [ya(1) - y0(1);yb(1) - y0(2);ya(2) - y0(3);yb(2) - y0(4)];
end

[Torsten]

I don't know what the size of lambda_0 is.

[KM]

lambda_0 is just a constant.

[Torsten]

Then you should remember what changes you made to the function because this one worked:

  dy1 = [y(3)...
         y(4)...
         y(3)./xi + (g^2 * (1+y(1)) * (1+(lambda^2*y(4)^2)) * sin(y(2))^2)...
         (1./(1+(lambda^2*(n*(y(1)+1)./xi).^2*sin(y(2))^2))) .* ( ((sin(y(2))*cos(y(2))*(n*(y(1)+1)./xi).^2) + (m^2*sin(y(2)))) - (y(4)./xi).*( ((lambda^2*(n*(y(1)+1)./xi)*sin(y(2))^2).*((n*(y(1)+1)./xi)+(2*xi.*(((xi*n*y(3))-(n*y(1))-n)./xi.^2)))) + 1 + (lambda^2*y(4)*xi.*(n*(y(1)+1)./xi).^2.*sin(y(2))*cos(y(2))) ) ) ];
  

Maybe because the ... are missing in the third line ?

And remember that in your new code, you use cot(x) which is Inf at all multiples of pi. This can easily lead to a singular Jacobian.