clear all

global a K

a=1;
K=10;

dt=0.28;
Tfinal=105;
Nt=ceil(Tfinal/dt);

xt1=zeros(Nt,1);xt4=xt1;xt3=xt1;tt=xt1;

x1=0.1;x3=x1;x4=x1;
xt1(1)=x1;xt4(1)=x4;xt3(1)=x3;tt(1)=0;

for ia=2:Nt;
    tt(ia)=(ia-1)*dt;
    %Euler
    k1=fdot1d(x1);
    x1=x1+dt*k1;  
    xt1(ia)=x1;
    
    %A3
    k1=fdot1d(x3);
    k2=fdot1d(x3+k1*dt/2);
    k3=fdot1d(x3+k2*3*dt/4);
    x3=x3+dt*(2*k1+3*k2+4*k3)/9;
    xt3(ia)=x3;
    
    %RK4
    k1=fdot1d(x4);
    k2=fdot1d(x4+k1*dt/2);
    k3=fdot1d(x4+k2*dt/2);
    k4=fdot1d(x4+k3*dt);
    x4=x4+dt*(k1+2*k2+2*k3+k4)/6;
    xt4(ia)=x4;
    
%     figure(32)
%     plot(tt(1:ia),xt(1:ia),'r',tt(1:ia),0*tt(1:ia),'k','linewidth',2)
%     xlabel('t');ylabel('x(t)')
%     drawnow
    

end;
     figure(32)
     plot(tt,xt1,'k',tt,xt3,'g',tt,xt4,'r',tt,0*tt,'k','linewidth',2)
     xlabel('t');ylabel('x(t)')
     %axis([0 1 0 11])
     drawnow
     %hold on

     figure(320)
     plot(tt,abs(xt1-xt4),'k',tt,abs(xt3-xt4),'r',tt,0*tt,'k','linewidth',2)
     xlabel('t');ylabel('x(t)')
     %axis([0 1 0 11])
     drawnow