function [t,x,epsilon] = RKAdaptive45(expression,int,x0,n,itmax,epsmin,epsmax,hmin,hmax)
% OR:
%function [t,x,epsilon] = RK45_Adaptive(int,x0,n,itmax,epsmin,epsmax,hmin,hmax)
% The ADAPTIVE Runge-Kutta-Fehlberg method of order 45 
%Input: Interval [t0,T], initial condition x0, 
%       n=the number of subintervals (determines the initial step size h)
%Output: the approximate solution at time T, epsilon = local truncation
%error
%Example Usage: [t,x,epsilon] = RKAdaptive45('2.0 + (x - t -1.0).^2',[1 1.5625],2,72,1.0e3, 1.0e-8, 1.0e-5,1.0e-6,1.0)
% OR;
%Example Usage: [t,x,epsilon] = RKAdaptive45([1 1.5625],2,72,1.0e3, 1.0e-8, 1.0e-5,1.0e-6,1.0)

%n = 72;
ta = int(1); tb = int(2);
format long
%n = round((int(2)-int(1))/h);
h = (int(2)-int(1))/n;
t = int(1);
x = x0;
%itmax = 1000
delta = .5e-5;
iflag = 1;
k = 0;

c20 = .25;  c21 = .25;
c30 = .375; c31 = .09375;  c32 = .28125;
c40 = 12./13.; c41 = 1932./2197.;
c42 = -7200./2197.; c43 = 7296./2197;
c51 = 439./216.; c52 = -8.;
c53 = 3680./513.; c54 = -845./4104.;
c60 = .5; c61 = -8./27.; c62 = 2;
c63 = -3544./2565.; c64 = 1859./4104.;
c65 = -0.275;
a1 = 25./216.; a2 = 0.; a3 = 1408./2565.;
a4 = 2197./4104.; a5 = -0.2;
b1 = 16./135.; b2 = 0.; b3 = 6656./12825.;
b4 = 28561./56430.; b5 = -0.18;
b6 = 2./55.;
    
while k <= itmax
    k = k+1;
    if abs(h) < hmin
        h = sign(h)*hmin;
    end
    if abs(h) > hmax
        h = sign(h)*hmax;
    end
    d = abs(tb - t);
    if d <= abs(h)
        iflag = 0
        if d <= delta*max(abs(tb), abs(t)) break
            h=sign(h)*d;
        end
    end

    x_save=x;
    t_save=t;

   
    % the begining of the RK4 procedure
    
    K1 = h*xdot(t,x);
    K2 = h*xdot(t+c20*h,x+c21*K1);
    K3 = h*xdot(t+c30*h,x+c31*K1 + c32*K2);
    K4 = h*xdot(t+c40*h,x+c41*K1+c42*K2+c43*K3);
    K5 = h*xdot(t+h,x+c51*K1+c52*K2+c53*K3+c54*K4);
    K6 = h*xdot(t+c60*h,x+c61*K1+c62*K2+c63*K3+c64*K4+c65*K5);    
    x4 = x + a1*K1 + a3*K3 + a4*K4 + a5*K5;
    x = x + b1*K1 + b3*K3 + b4*K4 + b5*K5 + b6*K6;
    t = t + h;
    epsilon = norm(x-x4);
    % the end of the RK4 procedure    

    if iflag < 1 break   %EXIT WHILE LOOP HERE
        %= 0 break   %EXIT WHILE LOOP HERE
    end
    if epsilon < epsmin
        h = 2*h
    end
    if epsilon > epsmax
        h = h/2;
        x = x_save;
        t = t_save;
        k = k-1;
    end
figure(1)    
plot(t,x,'.','LineWidth',1,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);
hold all
legend('RK45')

    
end %end of the while loop


function z = xdot(expression,t,x)
syms t x
z = subs(expression,x,t);
%function z = xdot(t,x)
%z = 2.0 + (x - t -1.0).^2;