%Problem 1 of page 178 of  Turner
% General solution to y' = x/y
clc
clear all
close all
yprime=inline('x/y','x','y');
%General solution
ygen = vectorize(inline('sqrt(x*x + b)','x','b'));
xlo=0; xhi=4;
xvec = linspace(xlo,xhi,100);
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.0;
Steps = [10 20 50 100 200];
symbols = ['*-r';'o-g';'<-b';'s-y';'x-m'];
hold on
for index = 1:5
    Nsteps = Steps(index);
    h=(xhi-xlo)/Nsteps;
    alfa = .5;
    XY=eulerMod(yprime,xlo,xhi,y0,Nsteps,alfa);
    %plot(xvec,yvec,'k',x1,y1,'*-r')
    plot(XY(:,1),XY(:,2),symbols(index,:));
    hsize(index)=h;
    error(index) = abs(XY(end,2) - ygen(4,9));
end;
plot(xvec,yvec,'k')
hold off
xlabel ('x')
ylabel ('y(x)')
title ('Euler Solution: Nsteps = 10, 20, 50, 100, and 200')
hsize
error
error./hsize
%legend(bvec,'Location','EastOutside')
