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

clear
% ************************
% User Inputs
stepsize= 0.01;
f0= 0.0;   % initial condition at xI [i.e. f(xI)]
xI= 0;  % boundary conditions
xF= pi/2;

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

F(1)= f0;  %initialize first value
k=1;  % counter
for i=xI:stepsize:xF
 if i == xI
  F(k)= f0;
 else
  F(k)= F(k-1) + stepsize* (F(k-1)^2 + 1);  % first calculate f at this step
  nextstep= (F(k) + stepsize*F(k))^2 + 1;   % use that to get it at the subsequent step
  % now use these two to effectively average what f is over this interval
  F(k)= F(k-1) + (stepsize/2) * (F(k-1)^2 + 1 + nextstep);
 end
 x(k)= i;  % 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= f^2 +1 using improved Euler method')

% plot analytic solution
sol= tan(x);
plot(x,sol, 'r.-')
legend('Numeric solution', 'True solution', 0); axis([0 pi/2 0 100])