% York U PHYS 2030  (09.15.14)
% use basic Euler method to solve the following differential equation
% f'(x)= x^2 + 1

clear
% ************************
% User Inputs
Fprime = @(x)(x.^2+1); % function describing the ODE
Fsol= @(x,IC)(1/3*x.^3 + x + IC); % solution (known; IC is the initial condition, f0)
stepsize= 0.5;
f0= 0.0;   % initial condition at xI [i.e. f(xI)]
xB= [0 2];  % interval to solve over (i.e., boundaries)

% ************************

F(1)= f0;  %initialize first value
k=1;  % counter

for j=xB(1):stepsize:xB(2)
 if (j == xB(1)),   
     F(k)= f0;  % handle initial condition
 else
     F(k)= F(k-1) + stepsize*(Fprime(j));    % apply Euler's method
 end
 x(k)= j;  % keep track of x for plotting
 k= k+ 1;  % increment counter
end

% visualize
figure(1); clf;
plot(x,F, 'o--','LineWidth',2);    grid on;    hold on
xlabel('x');    ylabel('f(x)'); title('Solution to df/dx= x^2 +1')

% now plot exact solution and derivative
plot(x,Fsol(x,f0), 'r.-','LineWidth',1);
plot(x,Fprime(x), 'kd-');
legend('f(x) (numeric solution)', 'f(x) (true solution)', 'df/dx (from ODE)','Location','NorthWest')