function nma_RK4(rate,stepSize,startX,endX,initialY)
%solve 1st order ODE using Runge-Kutta classical 4th order

%function nma_RK4(rate,stepSize,initialX,finalX)
%
% Function to solve 1st order ODE using Runge-Kutta classical 4th order
%
%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.
%
% EXAMPLE:
% rate='x^2'; stepSize=.1; startX=0;
% endX=1; initialY=1; nma_RK4(rate,stepSize,startX,endX,initialY)
% step	 x		 y(x) 	  end_step_y
% 1	  0.0000	  1.0000	  1.0003
% 2	  0.1000	  1.0003	  1.0027
% 3	  0.2000	  1.0027	  1.0090
% 4	  0.3000	  1.0090	  1.0213
% 5	  0.4000	  1.0213	  1.0417
% 6	  0.5000	  1.0417	  1.0720
% 7	  0.6000	  1.0720	  1.1143
% 8	  0.7000	  1.1143	  1.1707
% 9	  0.8000	  1.1707	  1.2430
% 10	  0.9000	  1.2430	  1.3333
% 11	  1.0000	  1.3333	  1.3333
%
% copyright: Nasser M. Abbasi

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

currentPoint = 1;
data = zeros(nPoints,4);

while(currentPoint < nPoints)
    
    startStepX =  x;
    startStepY =  y;
    
    k1 = yder(rate,startStepX           ,startStepY );
    k2 = yder(rate,startStepX+stepSize/2,startStepY+(1/2)*stepSize*k1 );
    k3 = yder(rate,startStepX+stepSize/2,startStepY+(1/2)*stepSize*k2 );
    k4 = yder(rate,startStepX+stepSize  ,startStepY+stepSize*k3 );
    
    nStep         = nStep+1;
    data(nStep,1) = nStep;
    data(nStep,2) = startStepX;
    data(nStep,3) = startStepY;
    y             = startStepY + stepSize*( (1/6) * (k1+ 2*k2+ 2*k3+ k4) );
    data(nStep,4) = y;
    currentPoint = currentPoint+1;
    x = x + stepSize;
end

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

analyze(data);

%my_test(data,stepSize);
end

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

%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%
function analyze(d)

[row,~]=size(d);
fprintf('step\t x\t\t y(x) \t  end_step_y\n');


for i=1:row
    fprintf('%d\t %7.4f\t %7.4f\t %7.4f\n',...
        d(i,1),d(i,2),d(i,3),d(i,4));
end

end

%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%
function my_test(d,stepSize)
plotAnalyticAndRKOnly(d,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,2);
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(d(:,2),d(:,3),'.:')

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','RK4 h=0.5');
ylabel('y');
title(sprintf('analytic solution of %s compared to RK4 (h=%.2f)',...
    y_analytic,stepSize));
end

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

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

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