%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Solves the IVP y' = sqrt(2*a^2*(y+e^-y) +c^2*cos(y)^2), y(0)=1
%on the time interval [0,T] using the explicit Euler method and 
%Heun's method.  Calls the function f, which defines the RHS of 
%the equation.
%
%Input Variables:   a, c, T, 
%                   nts: number of time steps, entered by the user
%
%Output Variables:  SE, SH: solution vectors for Euler and Heun resp.
%                   EE, EH: Richardson error estimate for Euler and Heun
%                   Efeval, Hfeval: # of func. eval. for Euler and Heun
%
%Variables Used in Computation: 
%                   h: step size
%                   t: time vector
%                   u, U: iterate values for Euler
%                   v, V: iterate values for Heun
%
%Printed Output: The value of h used in the computation, along 
%               with the number of function evaluations and the 
%               maximum error estimate for each method.
%Plotted Output: Numerical solutions and corresponding error 
%               estimates as functions of time.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear

% INITIAL DATA
a = 1; c = 3;
nts = input ('number of time steps: ');
T = 10; h = T/nts;
t = 0:h:T+2*h;
u = 1; v = u; U = u; V = v;
Efeval = 0; Hfeval = 0;

% INTEGRATION LOOP
for j = 1:2:nts+1

   %EULER'S METHOD
   SE((j+1)/2) = u;
   EE((j+1)/2) = u -U;
   for k = j:j+1
        u = u + h*f(t(k),u,a,c);
   end
   U = U + 2*h*f(t(j),U,a,c);
   Efeval = Efeval +3;
    
   %HEUN'S METHOD
   SH((j+1)/2) = v;
   EH((j+1)/2) = (v -V)/3;
   for k = j:j+1
        fv = f(t(k),v,a,c);
        v = v + (h/2)*(fv +f(t(k+1),v +h*fv,a,c));
   end
   fV = f(t(j),V,a,c);
   V = V + h*(fV +f(t(j+2),V +2*h*fV,a,c));
   Hfeval = Hfeval +6; 
end

% Print results 
fprintf('step size = %e\n', h);
fprintf('\n');
fprintf('%i function evaluations were required for Euler.\n',Efeval);
fprintf('%e is the estimated error in Euler`s approximation.\n', norm(EE,inf));
fprintf('\n');
fprintf('%i function evaluations were required for Heun.\n',Hfeval);
fprintf('%e is the estimated error in Heun`s approximation.\n', norm(EH,inf));

% Plot results
t = 0:2*h:T;
subplot(2,1,1)
plot(t,SE,'b.-',t,SH,'r')
xlabel('t'); ylabel('y(t)')
subplot(2,2,3)
plot(t,EE,'b.-',t,EH,'r')
xlabel('t'); ylabel('error'); legend('Euler','Heun','Location','NorthWest')
subplot(2,2,4)
plot(t,EH,'r')
xlabel('t'); ylabel('error'); legend('Heun','Location','NorthWest')