function [node,dt] = varRK(node,elem,mg,mp,M,gamma,ki,dt,acc,maxdt)
%[node,delt] = varRK(solvertype,node,elem,m,M,gamma,ki,dt,acc)
%
%  node: Nodal positions
%  elem: The element connectivities
%  delt: The time step used
%   m,M: Mobilities
% gamma: The surface tension
%    ki: Used to determine the amount of anisotropy
%    dt: The initial time step to use(input), the actual time step
%        used(output)
%   acc: The prescribed accuracy

done=0;

count=0;

onemore=0;

while done==0

    count=count+1;
    
    %Find the updated nodal values at dt
    [vnodal]=fea(node,elem,mg,mp,M,gamma,ki);
    k1=vnodal*dt;
    [vnodal]=fea(node+k1/2,elem,mg,mp,M,gamma,ki);
    k2=vnodal*dt;    
    [vnodal]=fea(node+k2/2,elem,mg,mp,M,gamma,ki);
    k3=vnodal*dt;    
    [vnodal]=fea(node+k3,elem,mg,mp,M,gamma,ki);
    k4=vnodal*dt;
    nodetest{1}=node+(1/6)*(k1+2*k2+2*k3+k4);

    %Find the updated nodal values at t+dt taking two steps, dt/2 each.
    nodetest{2}=node;
	for i=1:2
        [vnodal]=fea(nodetest{2},elem,mg,mp,M,gamma,ki);
        k1=vnodal*dt/2;
        [vnodal]=fea(nodetest{2}+k1/2,elem,mg,mp,M,gamma,ki);
        k2=vnodal*dt/2;
        [vnodal]=fea(nodetest{2}+k2/2,elem,mg,mp,M,gamma,ki);
        k3=vnodal*dt/2;
        [vnodal]=fea(nodetest{2}+k3,elem,mg,mp,M,gamma,ki);
        k4=vnodal*dt/2;
        nodetest{2}=nodetest{2}+(1/6)*(k1+2*k2+2*k3+k4);
	end
	
    %Find the absolute maximum error
	delta=max(abs(nodetest{2}-nodetest{1})); %The maximum errors in the x and y directions
    delta=max(delta); %The overall maximum error
    
    %If the maximum error is bewteen 10% of the required accuracy and the
    %accuracy, then this process is done. If not, find a new time step and
    %run this again.
    if onemore==0;
		if and(or(and(delta<=acc,delta>=0.1*acc),delta==0),dt<=maxdt) %just right
            node=nodetest{2};
            done=1;
        elseif and(and(count>25,delta<=acc),dt<=maxdt) %It might never find a value that will satisfy the above condition
            if and(dt>maxdt,isempty(maxdt)==0)
                dt=maxdt;
            else
                node=nodetest{2};
                done=1;
            end
        elseif and(and(dt>maxdt,~isempty(maxdt)),delta<=acc)
            %If it is accurate enough, but the maximum time step has been
            %surpased, then use the maximum time step.
            dt=maxdt; 
            onemore=1;
        else
            dt=dt*(acc/delta)^(0.2);
        end
    else
        done=1;
    end

end