How do you solve a boundary value problem of the given form ?

The differential equation I want to solve is d2y/dx2 = k*(y^n) , where k and n are constants and n<0. The boundary conditions are y=1@x=1 and dy/dx=0@x=0. I tried using the bvp4c solver for this case but it always showed " Error using bvp4c (line 251) : Unable to solve the collocation equations -- a singular Jacobian encountered." I want to know if there is any alternative way to solve this equation. Thanks.

[Walter Roberson]

Is n a negative integer, or is it floating point? If it is floating point, then which meaning of raising to a power do you want to use for the negative y(x) ?

[Torsten]

Please show us the MATLAB code you are using.

Best wishes

Torsten.

[Anshuman Goswami]

Here's the code as requested. And y(x) is always positive, since I am only solving for x between 0 and 1. 'n' can be a floating point as well as an integer. In the code below, I have taken n = -0.25 but it may be -1 or -2 as well.

---------------------------------------------------------------------

%main.m
xlow = 0; 
xhigh = 1; 
global phi ncall phi2 n;
n = -0.25;
ncall = 1; 
phi = [0.01 0.5 1.0 5.0 10.0] ; 
for ncall=1:5
solinit = bvpinit(linspace(xlow,xhigh,10),[0 0]);
sol = bvp4c(@bvp4ode,@bvp4bc,solinit);
xint = linspace(xlow,xhigh);
Sxint = deval(sol,xint);
if ncall==1
p = plot(xint,Sxint(1,:),'c');
else if ncall==2 
p = plot(xint,Sxint(1,:),'m');
else if ncall==3
p = plot(xint,Sxint(1,:),'r');
else if ncall==4 
p = plot(xint,Sxint(1,:),'g');
else if ncall==5 
p = plot(xint,Sxint(1,:),'b');
end
end
end
end
end
end

--------------------------------------------------------------------

%bvp4ode.m
function dydx = bvp4ode( x,y)
%BVP4ODE Summary of this function goes here
%   Detailed explanation goes here
global phi ncall n
dydx = [ y(2) phi(ncall)*phi(ncall)*(y(1)^n)];
end

---------------------------------------------------------------------

%bvp4bc.m
function res = bvp4bc( ya,yb )
%BVP4BC Summary of this function goes here
%   Detailed explanation goes here
res = [ ya(2) yb(1)-1 ] ; 
end

-------------------------------------------------------------------

Answer 1

[Walter Roberson]

Your lower bounds appear to be xlow = 0, which you are raising to a negative power. That is going to generate 1/0 which is a singularity. Your xlow needs to be at least eps(realmin)

Answer 2

[Torsten]

Replace

solinit = bvpinit(linspace(xlow,xhigh,10),[0 0]);

by

solinit = bvpinit(linspace(xlow,xhigh,10),[1 0]);

Best wishes

Torsten.