You did not initialize y_h(1) before you used it.
Your plots are very different scales: your analytic goes up to about 40 but your heun goes to about 10^10 .
This suggests that you have programmed your heun incorrectly.
format long g
n = 127;
x = linspace(0,2,n);
h = x(n)-x(n-1);
y_a =((x.^2+1/2.*x+1)).^2;
y_h(1) = 1;
for i=1:n-1
A(i) = x(i) + h;
B(i) = y_h(i) +h*(1+4.*x(i)).*sqrt(y_h(i));
L(i) = (1+4.*x(i)).*sqrt(y_h(i));
R(i) = (1+4.*A(i)).*sqrt(B(i));
y_h(i+1) = y_h(i) + 1/2.*(h.*L(i).*R(i));
end
plot(x,y_a,'--kd')
figure
plot(x,y_h,'--rx')
