%Problem 8 of page 186 of  Turner
% General solution to y' = x + y*y
clc
clear all
close all
yprime=inline('x/y','x','y');
%General solution
ygen = vectorize(inline('sqrt(x*x + b)','x','b'));
%yprime=inline('x+y*y','x','y');
%General solution
%ygen = vectorize(inline('(b*airy(1,-x)+airy(3,-x))/(b*airy(0,-x)+airy(2,-x))','x','b'));
xlo=0; xhi=1;
xvec = linspace(xlo,xhi,100);
sqrt3 = sqrt(3);
b = 5;
indexhi=10;
for index = 1:indexhi
    yvec(index,:)= ygen(xvec,b);
   % bvec(index,:) = num2str(b); 
    bvec{index} = num2str(b);
    b = b+1;
end
%numerical Solution
y0=3;
%Steps = [10 20 50 100 200];
symbols = ['*-r';' -g';' -b';'s-y';'x-m'];
hold on
for index = 1:3
    h= 10^(-index);
     Nsteps = (xhi-xlo)/h;
    XY=RK4(yprime,xlo,xhi,y0,Nsteps);
    %plot(xvec,yvec,'k',x1,y1,'*-r')
    plot(XY(:,1),XY(:,2),symbols(index,:));
    hsize(index)=h;
    error(index) = abs(XY(end,2) - ygen(xhi,9));
end;
plot(xvec,yvec,'k')
hold off
xlabel ('x')
ylabel ('y(x)')
title ('RK4 Solution:y'' = x/y: Nsteps = 10, 100, 1000')
format long
hsize
error
error./hsize
%legend(bvec,'Location','EastOutside')
