function nma_euler_heun(rate,stepSize,startX,endX,initialY,tol)
% Solve ODE using Euler-Heun (corrector-predictor method)

%INPUT
%   rate: A string that represents the derivative dy/dx. This
%   is the function f(x,y) in the ODE equation dy/dx= f(x,y)
%   This string MUST use the letters 'x' and 'y' for the
%   independent and the dependent variables respectively.
%   For an example, in the equation dy/dx = y*x, the rate
%   string is passed in as 'y*x'
%
%   stepSize: This is 'h', the step size over x
%
%   startX: The starting value of x
%   endX:   The ending value of x
%   initialY: Initial condition. The value of y at startX.
%   tol: In percentage. The value to stop the corrector-predictor
%        loop at each step.
%
% EXAMPLE
% rate='x^2'; stepSize=.1; startX=0; endX=1; initialY=1;
% nma_euler_heun(rate,stepSize,startX,endX,initialY,.001)

%
% copyright: Nasser M. Abbasi
% 2003

% This is protect again predictor-correctpr loop not converging
%to with epsilon.
MAX_STEPS_PER_ONE_ITERATION=1000;

y = initialY;
x = startX;
nStep   = 0;
nPoints = (endX-startX)/stepSize +1;

currentPoint = 1;

while(currentPoint < nPoints)
    
    converged  =  false;
    
    startStepX =  x;
    endStepX   =  x + stepSize;
    startStepY = y;
    
    rateStartOfStep = yder( rate, startStepX, startStepY );
    
    y_predict       = startStepY + ( stepSize * rateStartOfStep );
    
    thisIterStep    = 0;
    
    while(~converged)
        
        nStep         = nStep+1;
        data(nStep,1) = nStep;
        data(nStep,2) = startStepX;
        data(nStep,3) = startStepY;
        data(nStep,4) = y_predict;
        rateEndOfStep = yder( rate, endStepX, y_predict );
        
        averageRate   = ( rateStartOfStep + rateEndOfStep )/2;
        
        y_correct     = startStepY + ( stepSize * averageRate );
        
        data(nStep,5) = y_correct;
        
        relativeError = abs( (y_correct - y_predict)/y_correct )*100;
        
        data(nStep,6) = relativeError;
        data(nStep,7) = averageRate;
        
        y_predict     = y_correct;
        thisIterStep  = thisIterStep+1;
        
        if(relativeError<tol || thisIterStep>MAX_STEPS_PER_ONE_ITERATION)
            converged=true;
        end
    end
    
    % Update for next point
    x  =  x+stepSize;
    y  =  y_correct;
    
    currentPoint = currentPoint+1;
end

nStep         = nStep+1;
data(nStep,1) = nStep;
data(nStep,2) = x;
data(nStep,3) = y;
data(nStep,4) = y;
data(nStep,5) = y;

data(nStep,6) = 0;
data(nStep,7) = 0;

analyze(data,stepSize);


%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%
function dydx=yder(f,x,y)
dydx=eval(f);

%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%
function analyze(d,stepSize)

[row,~]=size(d);
fprintf('step\t x\t\t y(x) \t  predict_y(x+h) correct_y(x+h) \tepsilon(a)\n');

N=0;
for i=1:row
    if(i>1)
        if(d(i-1,2) == d(i,2))
            fprintf('%d\t                         %7.4f\t %7.4f\t\t %7.4f\t%7.4f\n',...
                d(i,1),d(i,4),d(i,5),d(i,6),d(i,7));
        else
            N=N+1;
            cleanData(N,1)=d(i-1,2);
            cleanData(N,2)=d(i-1,3);
            
            fprintf('%d\t %7.4f\t %7.4f\t %7.4f\t %7.4f\t\t %7.4f\t%7.4f\n',...
                d(i,1),d(i,2),d(i,3),d(i,4),d(i,5),d(i,6),d(i,7));
            
        end
    else
        fprintf('%d\t %7.4f\t %7.4f\t %7.4f\t %7.4f\t\t %7.4f\t%7.4f\n',...
            d(i,1),d(i,2),d(i,3),d(i,4),d(i,5),d(i,6),d(i,7));
    end
    
    if(i<row)
        if(d(i,2) ~= d(i+1,2))
            fprintf('\n');
        end
    end
end

N=N+1;
cleanData(N,1)=d(end,2);
cleanData(N,2)=d(end,3);

cleanData;

%my_test(cleanData,stepSize)

%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%
function my_test(cleanData,stepSize)
plotAnalyticAndHeunOnly(cleanData,stepSize);

figure;

x_5=[0 0.5 1 1.5 2];
y_5=[1 0.4 0.21 0.189 0.288225];
x_25=[0 0.25 0.5 0.75 1 1.25 1.5 1.75 2];
y_25=[1  0.7 0.5009375  0.38196484375     0.321089196777344 ...
    0.305034736938477  0.332678509973526  0.332678509973526 0.487581941179949];

y_analytic='exp(x*(x^2/3 - 1.2))';
ezplot(y_analytic,0,3);
hold on;
plot(x_5,y_5,'o')
plot(x_5,y_5,'-.')
plot(x_25,y_25,'x')
plot(x_25,y_25,'--')

plot(cleanData(:,1),cleanData(:,2),'.:')

xlim([-0.1 2.1]);
ylim([0 2.2])
legend('analytic','h=0.5 pts','h=0.5','h=0.25 pts','h=0.25','Heun h=0.25');
ylabel('y');
title(sprintf('analytic solution of %s compared to Euler and Heun(h=%.2f)',y_analytic,stepSize));


%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%
function plotAnalyticAndHeunOnly(cleanData,stepSize)
figure;

y_analytic='exp(x*(x^2/3 - 1.2))';
ezplot(y_analytic,0,4);
hold on;

plot(cleanData(:,1),cleanData(:,2),'.:')
xlim([-0.1 2.1]);
ylim([0 2.2])
legend('analytic','Heun h=0.25');
ylabel('y');
title(sprintf('analytic solution of %s compared to Heun(h=%.2f)',y_analytic,stepSize));