Unable to solve boundary condition problem

Hello, I've got an ODE with a singularity at x=0. I'm trying to solve for D2y but am having difficulty. Is there anyone who can help me with this problem? Thanks

syms y Dy D2y x  
k=5/8;
w=10;
q=0.5;
p1=1;
r=0.05;
e=2;
a=0.3;
t=0;
s1=2;
m2=0.5;
f=0.05;                    
o=1/14;
z=0.2;
vmax=100;
sigma=(1-k)*(e/(1-o*e));
xf=2200;
b=1-((r+z*f*(a*e/(1-a)*(1-o*e)+1)/f*(a*e/(1-a)*(1-o*e))));
m=(Dy*x*(sigma*m2+0.5*sigma*(sigma-1)*(s1^2))/y)+(0.5*D2y*(x^2)*(sigma^2)*(s1^2)/(2*y));
s=Dy*sigma*s1/y;
t1=(xf*(e/(e-1))* (1/a-1)^(e/(1-o*e)-1)*q^(1-e/(1-o*e))/w^(e/(1-o*e))*((z*f/(r+t*p1*s1-m2))+1)*(1/(p1*s1+r+(1-b)*f-m2))^((a*e/(1-a)*(1-o*e))+1));
ode=y-(x*(e/(e-1))* (1/a-1)^(e/(1-o*e)-1)*q^(1-e/(1-o*e))/w^(e/(1-o*e))*((z*f/(r+t*p1*s-m))+1)*(1/(p1*s+r+(1-b)*f-m))^((a*e/(1-a)*(1-o*e))+1));
D2ynum = solve(ode==0,D2y);
f = matlabFunction(D2ynum);
odefcn = @(x,y)[y(2);f(y(2),x,y(1))];
bcfcn = @(ya,yb)[ya(1);yb(1)-t1];
xmesh = linspace(0.001,xf,10);
solinit = bvpinit(xmesh, [1 1]);
sol = bvp4c(odefcn,bcfcn,solinit);

Answer 1

[Torsten]

"ode == 0" cannot be explicitly solved for D2y. But given y and Dy, you can try to use Newton's method to supply D2y.

Define a function (not a function handle) "odefcn" in which you call "fsolve" or "fzero" to solve ode == 0 for D2y given numerical values for y and Dy and return [y(2);result of ode == 0] to "bvp4c".

[Della]

I tried to solve ode==0 for D2y given numberical values for y and Dy but still can't derive D2y.

[Torsten]

Show us your attempt.

[Della]

syms y Dy D2y x  
k=5/8;
w=10;
q=0.5;
p1=1;
r=0.05;
e=2;
a=0.3;
t=0;
s1=2;
m2=0.5;
f=0.05;                    
o=1/14;
z=0.2;
vmax=100;
sigma=(1-k)*(e/(1-o*e));
xf=2200;
b=1-((r+z*f*(a*e/(1-a)*(1-o*e)+1)/f*(a*e/(1-a)*(1-o*e))));
m=(0.5*x*(sigma*m2+0.5*sigma*(sigma-1)*(s1^2)))+(0.5*D2y*(x^2)*(sigma^2)*(s1^2)/2);
s=0.5*sigma*s1;
t1=(xf*(e/(e-1))* (1/a-1)^(e/(1-o*e)-1)*q^(1-e/(1-o*e))/w^(e/(1-o*e))*((z*f/(r+t*p1*s1-m2))+1)*(1/(p1*s1+r+(1-b)*f-m2))^((a*e/(1-a)*(1-o*e))+1));
ode=1-(x*(e/(e-1))* (1/a-1)^(e/(1-o*e)-1)*q^(1-e/(1-o*e))/w^(e/(1-o*e))*((z*f/(r+t*p1*s-m))+1)*(1/(p1*s+r+(1-b)*f-m))^((a*e/(1-a)*(1-o*e))+1));
D2ynum = solve(ode==0,D2y);

[Torsten]

I plotted "ode" as a function of D2y with your initial values for x (e.g. x=0.001 in the below code), y(1) = 1 and y(2) = 1. But it does not have a zero (at least in the range -1500 <= D2y <= 1500). And "fzero" also gives up solving for D2y.

syms y Dy D2y x  
k=5/8;
w=10;
q=0.5;
p1=1;
r=0.05;
e=2;
a=0.3;
t=0;
s1=2;
m2=0.5;
f=0.05;                    
o=1/14;
z=0.2;
vmax=100;
sigma=(1-k)*(e/(1-o*e));
xf=2200;
b=1-((r+z*f*(a*e/(1-a)*(1-o*e)+1)/f*(a*e/(1-a)*(1-o*e))));
m=(Dy*x*(sigma*m2+0.5*sigma*(sigma-1)*(s1^2))/y)+(0.5*D2y*(x^2)*(sigma^2)*(s1^2)/(2*y));
s=Dy*sigma*s1/y;
t1=(xf*(e/(e-1))* (1/a-1)^(e/(1-o*e)-1)*q^(1-e/(1-o*e))/w^(e/(1-o*e))*((z*f/(r+t*p1*s1-m2))+1)*(1/(p1*s1+r+(1-b)*f-m2))^((a*e/(1-a)*(1-o*e))+1));
ode=y-(x*(e/(e-1))* (1/a-1)^(e/(1-o*e)-1)*q^(1-e/(1-o*e))/w^(e/(1-o*e))*((z*f/(r+t*p1*s-m))+1)*(1/(p1*s+r+(1-b)*f-m))^((a*e/(1-a)*(1-o*e))+1));
odenum = matlabFunction(ode)

odenum = function_handle with value:
    @(D2y,Dy,x,y)y+x.*(1.0./((Dy.*(7.0./4.0))./y-(D2y.*x.^2.*(4.9e+1./6.4e+1))./y-(Dy.*x.*(7.0./3.2e+1))./y+6.524468971261975e-2)).^(8.5e+1./4.9e+1).*(1.0./((D2y.*x.^2.*7.65625e+1)./y+(Dy.*x.*(1.75e+2./8.0))./y-5.0)-1.0).*7.239452177846195e-2

D2y = -1500:0.1:1500;
plot(D2y,odenum(D2y,1,0.001,1))

fun = @(D2y)odenum(D2y,1,0.001,1);
sol = fzero(fun,3)

Exiting fzero: aborting search for an interval containing a sign change
    because complex function value encountered during search.
(Function value at 2.84719e+06 is 0.99972+0.00030665i.)
Check function or try again with a different starting value.

sol = NaN