%Problem 1 of page 185 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=1;
xvec = linspace(xlo,xhi,100);
b = 0;
indexhi=15;
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 = [1 2 4];
symbols = ['*-r';'o-g';'<-b'];
hold on
for index = 1:3
    Nsteps = Steps(index);
    h=(xhi-xlo)/Nsteps;
    alfa=0.5;
    XY=eulerMod(yprime,xlo,xhi,y0,Nsteps,alfa);
    xeul = XY(:,1);
    yeul = XY(:,2);
    %plot(xvec,yvec,'k',x1,y1,'*-r')
    plot(xeul,yeul,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 ('Euler Solution: Nsteps = 1, 2, 4')
%legend(bvec,'Location','EastOutside')
hsize
error
error./hsize

