
function [] = prog1main()

dd=[10^(-5),10^(-8),10^(-10)];
for i=1:3
t=0;
T=17.0652166;
x=[0.994,0,0,-2.001585106];
tol=dd(i);
h=10^(-6);
   while t<T 
       
       b1=rkk(h,t,x);
       b2=rkkk(h,t,x);
       EST=nnorm(b1,b2);
       h=steuer(h,tol,EST);
        if (EST<=tol)
         x=rk(h,t,x);
         t=t+h;
         h=min(h,T-t);
        end    
   end
   
   disp(x)
  
end
    
function r=rk(tao,t,x)
 c=[0,1/5,3/10,4/5,8/9,1,1];
 b1=[35/384,0,500/1113,125/192,-2187/6784,11/84,0];
 b2=[5179/57600,0,7571/16695,393/640,-92097/339200,187/2100,1/40];
 A=[0,0,0,0,0,0,0;1/5,0,0,0,0,0,0;3/40,9/40,0,0,0,0,0;44/45,-56/15,32/9,0,0,0,0;19372/6561,-25360/2187,64448/6561,-212/729,0,0,0;9017/3168,-355/33,46732/5247,49/176,-5103/18656,0,0;35/384,0,500/1113,125/192,-2187/6784,11/84,0];
 
 k(1,1:4)=f(t,x);
 bk=b2(1)*k(1,1:4);
 for i=2:7
    he=[0,0,0,0]; 
     for j=1:(i-1)
         for ii=1:4
          he(ii)=he(ii)+A(i,j)*k(j,ii);
         end
     end    
    k(i,1:4)=f(t+c(i)*tao,tao*he);
    bk=bk+b2(i)*k(i,1:4);
 end

 r=x+tao*bk; 
 
function y=f(t,x)
  u=0.012277471;
  du=1-u;
  N1=((x(1)+u)^2+x(2)^2)^(3/2);
  N2=((x(1)-du)^2+x(2)^2)^(3/2);
  y=[x(3),x(4),x(1)+2*x(4)-du*((x(1)+u)/N1)-u*((x(1)-du)/N2),x(2)-2*x(3)-du*(x(2)/N1)-u*(x(2)/N2)];
  
 function rr=rkk(tao,t,x)
 c=[0,1/5,3/10,4/5,8/9,1,1];
 b1=[35/384,0,500/1113,125/192,-2187/6784,11/84,0];
 b2=[5179/57600,0,7571/16695,393/640,-92097/339200,187/2100,1/40];
 A=[0,0,0,0,0,0,0;1/5,0,0,0,0,0,0;3/40,9/40,0,0,0,0,0;44/45,-56/15,32/9,0,0,0,0;19372/6561,-25360/2187,64448/6561,-212/729,0,0,0;9017/3168,-355/33,46732/5247,49/176,-5103/18656,0,0;35/384,0,500/1113,125/192,-2187/6784,11/84,0];
 
 k(1,1:4)=f(t,x);
 bk=b1(1)*k(1,1:4);
 for i=2:7
    he=[0,0,0,0]; 
     for j=1:(i-1)
         for ii=1:4
          he(ii)=he(ii)+A(i,j)*k(j,ii);
         end
     end    
    k(i,1:4)=f(t+c(i)*tao,tao*he);
    bk=bk+b1(i)*k(i,1:4);
 end

 rr=tao*bk; 
 
 function rrr=rkkk(tao,t,x)
 c=[0,1/5,3/10,4/5,8/9,1,1];
 b1=[35/384,0,500/1113,125/192,-2187/6784,11/84,0];
 b2=[5179/57600,0,7571/16695,393/640,-92097/339200,187/2100,1/40];
 A=[0,0,0,0,0,0,0;1/5,0,0,0,0,0,0;3/40,9/40,0,0,0,0,0;44/45,-56/15,32/9,0,0,0,0;19372/6561,-25360/2187,64448/6561,-212/729,0,0,0;9017/3168,-355/33,46732/5247,49/176,-5103/18656,0,0;35/384,0,500/1113,125/192,-2187/6784,11/84,0];
 
 k(1,1:4)=f(t,x);
 bk=b2(1)*k(1,1:4);
 for i=2:7
    he=[0,0,0,0]; 
     for j=1:(i-1)
         for ii=1:4
          he(ii)=he(ii)+A(i,j)*k(j,ii);
         end
     end    
    k(i,1:4)=f(t+c(i)*tao,tao*he);
    bk=bk+b2(i)*k(i,1:4);
 end

 rrr=tao*bk; 
  
  function taonew=steuer(tao,tol,EST)
   
   taonew=tao*min(2.0,max(0.5,0.9*(tol/EST)^(1/(6+1))));
   
  function a=nnorm(bb1,bb2)
   a=0;
   for i=1:4
    a=a+(bb1(i)-bb2(i))^2;   
   end
   a=a^(1/2);