function [x, y] = eulerSolve(xstart, ystart, xend, h)
% eulerSolve: return data points using Euler's method to solve y' = xy^2 
%   xstart, ystart determine the initial condition
%   xend sets the end point and h is the step size
%   returns 2 column vectors, estimates of x and y(x) at n points
 
    % the "round" function gives the closest integer to some real number 
    n = round((xend-xstart+eps)/h);
    % adjust h so that the range is exactly n*h
    % (to ensure that we have exactly x(n+1) = xend
    hTrue = (xend-xstart)/n;

    x(1) = xstart;
    y(1) = ystart;
    
    for i=1:n
        yprime = x(i)*y(i)^2;
        y(i+1) = y(i) + hTrue*yprime;
        x(i+1) = x(i) + hTrue;
    end
return
