Singular problem for bvp4c

Hi All,

I encounter some singular problems regarding bvp4c, please see my code for my equations and implementation. I know that when x = -1 and 1, my derivative goes to infinity. And it seems like that MATLAB complaints about this. I wonder is there any workaround for this problem? I tried on ode113, and instead of evaluating my initial point at x = -1, I just do it at x = -1 + 1e-7 to avoid this singular points. And it is solvable. I would like to use bvp4c because I have other equations and boundary condition to coupled together. I wonder is there similar strategy for bvp4c? Here are my code

solinit = bvpinit(linspace(-1,1,50),0);
sol = bvp4c(@height_ode,@height_bc,solinit);
y = deval(sol,linspace(-1,1,50));
plot(linspace(-1,1,50),y(1,:))

and defintion for height_ode and height_bc are following:

function dhdx = height_ode(x,h)
w = 0.7;
q = 1/2*w*(1-x^2);
dqdx = -w*x;
kappa = 1.469;
dhdx = (560*h(1)^0.5*q - 64*h(1).^4 + 6*w.^(7/8)*(72*kappa - 77)*h(1)*q*dqdx)/(3*w.^(7/8)*(96*kappa - 77)*q.^2);
end

function res = height_bc(yleft, yright)
res = [yleft];

[Torsten]

And for bvp4c, the solution to integrate between [-1+eps:1-eps] does not work ?

Best wishes

Torsten.

[Zhekai Deng]

Thank you for your reply. I tried this approach, and it seems this strategy does not work, the height_ode are the same. To me, the only difference is that in ode113, the bcs is treated as an initial condition. In bvp4c, it is treated as a boundary condition on the left. See following code:

%%ode113 works
eps = 1e-6;
options = odeset('RelTol',1e-10);
xspan = [-1+eps 1-eps];
y0 =eps;
[x, h_head] = ode113(@height_ode,xspan,y0,options);
figure
plot(x,h_head,'-o')

%%bvp4c does not work
solinit = bvpinit(linspace(-1+eps,1-eps,3),eps);
sol = bvp4c(@height_ode,@height_bc,solinit);
y = deval(sol,linspace(-1+eps,1-eps,50));
plot(linspace(-1+eps,1-eps,50),y(1,:))

[Zhekai Deng]

Also, I noticed the initial condition I set for ode113 is not 0, but eps. However, even if I change height_bc as following, it still shows a singular problem. Any help would be greatly appreciated! Thank you

function res = height_bc(yleft, yright)
eps = 1e-5;
res = [yleft(1) + eps]; % res = [yleft(1) - eps]; also does not work 
end