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