2nd order diff equation for ode23/45?

I have spent hours trying to figure out why I'm not getting the correct answer.

%xy''-y'-8*x^3*y^3=0
%y''= y'/x + 8y^3x^2
% y'= a(1)
% y=a(2)
%y''(x) = a(1)/x + 8*(a(2)^3)*(a(1)^2)

y2p = @(x,a) [a(2); a(1)/x + 8*(a(2)^3)*(a(1)^2)];
tspan = [1 4];
y0 = [0.5 -0.5]

[y,yPrime]=ode23(y2p,tspan,y0)

%plot(y,yPrime)

[y2,y2Prime]=ode45(y2p,tspan,y0)

%plot(y2,y2Prime)

%exact solution
1/(1+1^2) %.5 at x=1
1/(1+4^2) %.0588 at x=4

Answer 1

[Davide Masiello]

The way you coded the ODE system is incorrect. See below

clear,clc

% Numerical solution
tspan       = [1 4];
y0          = [0.5 -0.5];
[x23,y23]   = ode23(@odefun,tspan,y0);
[x45,y45]   = ode45(@odefun,tspan,y0);

% Exact solution
x = linspace(1,4,100);
y = 1./(1+x.^2);

% Plot
plot(x,y,'k',x23,y23(:,1),'b',x45,y45(:,1),'r')
legend('exact','ode23','ode45')

% Function
function dydx = odefun(x,y) 
dydx(1,1) = y(2); 
dydx(2,1) = y(2)/x+8*y(1)^3*x^2;
end