Numerical solving second order differential equation

I would like to solve a second order differential equation,please see the attachment. I used the following code:

b = 0.19;
a = 0.072;
beta = 0.72;
Pr = 0.72;
c = beta/Pr*a^(4/3);

xspan = [0.01 100];
g40 = [1 0];
opts = odeset('RelTol',1e-4,'AbsTol',1e-6);
[x4,g4] = ode45(@(x4,g4) model4ode(x4,g4,b,c), xspan, g40, opts);

plot(log10(x4),g4)

function dg4dx4 = model4ode(x4,g4,b,c)
dg4dx4 = zeros(2,1);
dg4dx4(1) = g4(2);
dg4dx4(2) = (2./(3*x4)).*(4+x4.^(2*b).*(x4.^(2*b)+1).^(-1)).*g4(2)+...
    ((22/3)*c*x4.^(-2/3).*(x4.^(2*b)+1).^(1/(3*b))-(22/9)*x4.^(2*b-2).*(x4.^(2*b)+1).^(-1)).*g4(1);
end

But the result is not true. Could someone help me?

Answer 1

[Torsten]

b = 0.19;
a = 0.072;
beta = 0.72;
Pr = 0.72;
c = beta/Pr*a^(4/3);

xspan = [0.01 100];
g40 = [(xspan(1)^(2*b)+1)^(-1/(3*b))*(-11/3) 22/3*c*xspan(1)^(1/3)];
opts = odeset('RelTol',1e-4,'AbsTol',1e-6);
[x4,g4] = ode45(@(x4,g4) model4ode(x4,g4,b,c), xspan, g40, opts);
f4(:,1)=3/22*1/c*x4(:).^(-1/3).*g4(:,2);
f4(:,2)=(g4(:,1).*(x4(:).^(2*b)+1).^(1/(3*b))+0.5*1/c*x4(:).^(-1/3).*g4(:,2))./x4(:);
plot(log10(x4),f4)

function dg4dx4 = model4ode(x4,g4,b,c)
dg4dx4 = zeros(2,1);
dg4dx4(1) = g4(2);
dg4dx4(2) = (22/3*c*x4^(1/3)*(x4^(2*b)+1)^(1/(3*b))*g4(1)+g4(2)*(11+x4)/3)/x;
end

Best wishes

Torsten.

[Xiang Yi]

Thank you Torsten. The code works well, but the results are still not identical to the literature.

there is a bump, but we cannot obtain. Could you help me to find where is not correct?

[Star Strider]

What differential equation solver does ‘the literature’ use?

I would trust the MATLAB solver results.

[Torsten]

Could you include "the literature" as a pdf or a link ?

Best wishes

Torsten.

[Xiang Yi]

I'd like to solve the Eq.（20） in this paper. The parameters settings is the same as the Section 4, 1th paragraph. But we cannot obtain the results as Fig 1. and Fig. 2. You may note that there is a bump in the figure. x>0, and when x is very small and close to zero, f(x)=1; when x is very large and approaching positive infinity, f(x)=0.