function [tout, xout, nout] = VarStep(FuncF, tstep, x0, h0, tol)
%   VarStep	Integrates a system of ordinary differential equations using
%	a variable step variant of Mid-Point Euler.
%	[t,y,n] = VarStep('yprime', tspan, y0,h0,tol) integrates the system
%	of ordinary differential equations described by the M-file
%	yprime.m or inline function yprime over the interval tspan = [t0,tfinal] and using initial
%	conditions y0, h0, and tolerance tol.
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution vector.
%         yprime - Returned derivative vector; yprime(i) = dy(i)/dt.
% tspan = [t0, tfinal], where t0 is the initial value of t, and tfinal is
%         the final value of t.
% y0    - Initial value vector.
% h0    - Initial step size to try.
% tol   - Tolerance
%
% OUTPUT:
% t  - Returned integration time points (column-vector).
% y  - Returned solution, one solution row-vector per tout-value.
% n  - Returned number of evaluations of FuncF neccessary to compute the solution

% Initialization

t0=tstep(1);
tf=tstep(2);
i=1;

t=t0;
x=x0;
h = h0;

%Initialization of output

nout=0;
tout(i,1)=t;
xout(i,1)=x;


while(t<tf-10^-14)
    % Compute two approximations
    f1=feval(FuncF,t,x);
    f2=feval(FuncF,t+h/2,x+f1*h/2);
    A1=x+f1*h;
    A2=x+f1*(h/2)+f2*(h/2);
    
    nout=nout+2;
    
    % Estimate the Error Per Unit Step
    
    r=abs(A1-A2)/h;
   
    
    if(r>tol)
   % Decrease Step Size if neccessary
        h=min(tol/r *h,tf-t);
    else
    % Compute output result and iterate if satisfactory
        i=i+1;
        t=t+h;
        x=2*A2-A1;
        tout(i,1)=t;
        xout(i,1)=x;
        % Increase step size
        if(r > 0)
             h=min(tol/r*h,tf-t);
        else
             h=tf-t;
        end;

    end;

end;