Singular Jacobian Encountered Error- bvp4c

function sol = singular_flameballT1

format long e clc clear all global beta global gamma global Le global sl global S len1 = 0.0; % solving as a singular problem based on bvpset reference from matlab site len2 = 120;

beta = 10; gamma = 5; % gamma = q Le = 1.0; %DO NOT MODIFY THIS VALUE,CHANGE BELOW 4 lines sl = 0.977865;

x = len1:0.01:len2;

% generating the constant matrix for solving as a singular problem S = [0 0 0 -2.0];

options = bvpset('NMax',1200001,'SingularTerm',S,'FJacobian',@fjacob,'BCJacobian',@bcjacob,'Vectorized','on');

%load radLe1beta10len120step.mat;

solinit=bvpinit(x,@mat4init); sol=bvp4c(@eigen_value_sol,@bc,solinit,options); fprintf('First sol done!\n');

output(sol,len1,len2);

end

function [] = output(sol,len1,len2) global gamma global beta global Le global sl

figure(1) plot(sol.x,sol.y(1,:),'.');

hold on

Temp(1,:) = sol.y(1,:); Y(1,:) = 1-sol.y(1,:); plot(sol.x,Y(1,:),'r'); % calculating omega [~,b]=size(Temp); hold on m=zeros(1,b); for i=1:b m(1,i)=((1+gamma)/(1+(gamma*Temp(1,i)/Le)))*Y(1,i)*((beta^2)/(2*(sl^2)*Le))*exp((beta*(Temp(1,i)/Le-1))/((1+gamma*Temp(1,i)/Le)/(gamma+1))); end

fprintf('The omega at 0 is %13.9f.\n',m(1,1));

hold on %figure(3) plot(sol.x,m,'g'); % axis([len1 len2 -0.2 2]) legend('T','Y','w','location','NorthEastOutside'); xlabel('x') ylabel('solution') hold on

end

function f = eigen_value_sol(x,y)

global beta global gamma global Le global sl f=[ y(2,:) -((1+gamma)./(1+((gamma/Le)*y(1,:)))).*(1-y(1,:)).*(((beta^2)/(2*(sl^2)))*exp((beta*(y(1,:)/Le-1))./((1+(gamma/Le)*y(1,:))/(gamma+1)))) ]; end

function fjac = fjacob(x,y) global gamma global beta global sl global Le t1 = -((beta^2)/(2*(sl^2)))*(1+gamma); t3 = exp((beta*(y(1)/Le-1))/((1+gamma*y(1)/Le)/(gamma+1))); t2 = (1-y(1))/(1+(gamma*y(1)/Le)); dev1 = t1*t3*(1/(1+(gamma*y(1)/Le))^2)*(-1-gamma/Le) + t1*t2*t3*beta*(1+gamma)*(1/(1+(gamma*y(1)/Le))^2)*(1/Le)*(1+gamma);

fjac = [0 1 dev1 0];

end

function [bcajac,bcbjac] = bcjacob(x,y) % global gamma % global beta % global sl % global Le

bcajac = [0 1 0 0];

bcbjac = [0 0 1 0];

end

%boundary conditions function res=bc(ya,yb) %global Le %Dont change res=[ya(2) yb(1)];

end

%starting solution function yinit = mat4init(x) global Le c=4.0; a=2.0; % this is where the intersection happpens

% Le scaling and without derivatives yinit = [((1+erf(c*(a-x)))/2) 0];

% yinit = [((1+erf(c*(a-x)))/2) % ((-c/(pi^(1/2))*exp(-(c*(a-x))^2))) % ]; end

Any suggestions on why the singular jacobian error might occur in above code? Thanks in advance!