
clc % clear screen
clear % clear all variables

% this is similar to eulerExample but it shows that the approximate 
% method does not capture the divergence of y(x) at x=sqrt(C)
% ... it gives a finite estimate Y even as y(x) diverges
% ... as x->sqrt(C) the absolute error diverges, for all h.
%     useful to think about order of limits of h->0 and x->sqrt(C)...
xMin = 0;
xMax = 3;
yInitial = 1;
Cval = 2/yInitial + xMin*xMin;

% step size
h=0.1;
[xSol,ySol] = eulerSolve(xMin,yInitial,xMax,h);

% compute exact solution (as long as it exists)
nExact = 100;
for i=1:nExact
    xExact(i) = xMin + i*(xMax-xMin)/nExact;
    if ( xExact(i) < sqrt(Cval) )
        yExact(i) = 2/( Cval - xExact(i)*xExact(i) );
    else
        yExact(i) = NaN;  % "not a number", soln does not exist...
    end
end

% smaller step size
h=0.02;
[xSol2,ySol2] = eulerSolve(xMin,yInitial,xMax,h);

hold off 
plot( xExact,yExact, 'k-', 'Linewidth', 3 )
xlabel('x')
ylabel('y(x)')
xlim([ xMin xMax]);
% put a fixed upper limit (if reqd)
ylim([ 0 20 ])
set(gca,'FontSize',30)  % need to do this again for the new graph
hold on %% add the numerical solutions
plot( xSol,ySol, 'bo-' )
plot( xSol2,ySol2, 'mo-' )
legend('exact', 'h=0.1', 'h=0.02', 'Location','NorthWest')
hold off

% this command will print to a file, if req'd
% print -dpdf 'eulerFig2.pdf' -bestfit