Unable to solve the collocation equations -- a singular Jacobian encountered.

function Sc_1

% Define constants

phi = 0.02;

R_s = 1738;

R_f = 1053;

S_s = 230000000;

S_f = 0.18;

Cp_s = 1046.7;

Cp_f = 3594;

K_s = 156;

K_f = 0.492;

We = 0.3;

Ha = 0.3;

A = 0.1;

Pr = 4;

Q_star = 0.1;

Ec = 0.1;

s_1 = 0.1;

Sc = 0.3;

s_2 = 0.1;

Kr = 0.2;

s_3 =0.1;

Lb = 0.3;

Pe = 0.1;

delta_1 = 0.1;

k_prime = 0.2;

k_2prime = 0.1;

M_0 = 0.6;

alpha = pi/2;

n = -0.803;

B_2 = (1 - phi)^-2.5;

B_1 = (1-phi)+phi*(R_s/R_f);

B_4 = ((1-phi)+phi*((R_s*Cp_s)/(R_f*Cp_f)));

B_3 = ((S_s+2*S_f)-2*phi*(S_f-S_s))/((S_s+2*S_f)+phi*(S_f-S_s));

B_5 = ((K_s+2*K_f)-2*phi*(K_f-K_s))/((K_s+2*K_f)+phi*(K_f-K_s));

% Solve the BVP

% Create an options structure with specified tolerances and a Jacobian function

x = linspace(0, 1, 10);

options = bvpset('RelTol',1e-6,'AbsTol',1e-6);

solinit = bvpinit(x, [1 1 1 0 0 0 0 0 0]);

sol = bvp4c(@bvpexam2, @bcfun, solinit, options);

x_vals = sol.x;

y_vals = sol.y;

% Plot the solutions

figure(1);

plot(x_vals, y_vals(2,:), 'LineWidth', 1.3);

hold on; % Keep the plot for next iterations

figure(2);

plot(x_vals, y_vals(4,:), 'LineWidth', 1.3);

hold on; % Keep the plot for next iterations

figure(3);

plot(x_vals, y_vals(6,:), 'LineWidth', 1.3);

hold on; % Keep the plot for next iterations

figure(4);

plot(x_vals, y_vals(8,:), 'LineWidth', 1.3);

hold on; % Keep the plot for next iterations

% Boundary and ODE functions

function res = bcfun(ya, yb)

res = [Pr*ya(1) + (B_5/B_1)*M_0*ya(5)-0;

ya(2) - 1;

ya(4)+k_prime*ya(5)-0;

ya(6)+k_2prime*ya(7)-0;

ya(8)-0;

yb(2) - A;

yb(4) - (1 - s_1);

yb(6) - (1 - s_2);

yb(8) - (1 - s_3)];

end

function ysol = bvpexam2(x, y)

yy1 = (B_1*(y(2)^2 - y(1) * y(2)) + B_3*sin(alpha)^2*Ha*(y(2)-A) - A^2)/(B_2 + (3*(n- 1)/2)*B_2*We*y(3)*y(3));

yy2 = ((Pr*B_4)*(s_1*y(2) + y(2)*y(4) - y(1)*y(5)) ...

- B_2*Pr*Ec*((y(3))^2)*((1+(3*(n-1)/2))*We*(y(3))^2) + Q_star*(y(4)+s_1) ...

- B_3*Ha*Ec*Pr*((sin(alpha))^2)*((y(2) - A)^2))/B_5;

yy3 = Sc*(y(2)*y(6) + s_2*y(2) - y(1)*y(7)) - Kr*(y(6)+ 1/s_2);

yy4 = Lb*(y(8)+s_3) - Lb*y(9) + Pe*(y(9)*y(7) + (y(8)+delta_1+s_3));

ysol= [y(2); y(3); yy1; y(5); yy2; y(7); yy3; y(8); yy4];

end

end

[Torsten]

How does equation (9) in "Equations.png" come into play ?

[Farooq Aamir]

yy2 = ((Pr*B_4)*(s_1*y(2) + y(2)*y(4) - y(1)*y(5)) ...

- B_2*Pr*Ec*((y(3))^2)*((1+(3*(n-1)/2))*We*(y(3))^2) + Q_star*(y(4)+s_1) ...

- B_3*Ha*Ec*Pr*((sin(alpha))^2)*((y(2) - A)^2))/B_5;

This my equation number.

[Torsten]

Ah, I didn't see that it's the continuation of the preceeding line.

[Farooq Aamir]

@Torsten, Is it solovable?

[Torsten]

Is χ the same as N ?

And it should be

- B_2*Pr*Ec*((y(3))^2)*((1+(3*(n-1)/2))*We^2*(y(3))^2) + Q_star*(y(4)+s_1) ...

instead of

- B_2*Pr*Ec*((y(3))^2)*((1+(3*(n-1)/2))*We*(y(3))^2) + Q_star*(y(4)+s_1) ...

I guess.

There are no more differences I could find between code and mathematical description.

[Farooq Aamir]

@Torsten Still I have the same issue, Yes they are same.

[Torsten]

There can be many reasons that your code fails with this error message.

Mostly, it's due to programming errors, wrong parameter values or bad initial guesses for the functions.