clear all

global a e

a=0;e=0;
x0=0.2;

dt=0.1;

T=10;

Nt=ceil(T/dt);

x1=zeros(Nt,1);
tt=x1;

x1(1)=x0;
tt(1)=0;
for it=2:Nt
    %Euler
    k1=fdot1d(x1(it-1),tt(it-1));
    x1(it)=x1(it-1)+k1*dt;
    tt(it)=tt(it-1)+dt;
end;
   xan=x0./(x0+(1-x0)*exp(-tt));
   figure(1);hold on
   plot(tt,x1,'r',tt,exp(0*tt),'--k',tt,xan,'+b','linewidth',2)
   xlabel('t');ylabel('x');title('Euler')
   drawnow
   pause
   
x2=0*x1;
x2(1)=x0;
tt(1)=0;
for it=2:Nt
    %Feynman
    k1=fdot1d(x2(it-1),tt(it-1));
    k2=fdot1d(x2(it-1)+k1*dt,tt(it-1)+dt);
    x2(it)=x2(it-1)+(dt/2)*(k1+k2);
    tt(it)=tt(it-1)+dt;
end;
figure(1)
   plot(tt,x2,'g','linewidth',2)
   xlabel('t');ylabel('x');title('Euler,Feynman')
   pause
   
   
x3=0*x1;
x3(1)=x0;
tt(1)=0;

for it=2:Nt
    %Runge-Kutta
    k1=fdot1d(x3(it-1),tt(it-1));
    k2=fdot1d(x3(it-1)+k1*dt/2,tt(it-1)+dt/2);
    k3=fdot1d(x3(it-1)+k2*dt/2,tt(it-1)+dt/2);
    k4=fdot1d(x3(it-1)+k3*dt,tt(it-1)+dt);
    
    x3(it)=x3(it-1)+(dt/6)*(k1+2*k2+2*k3+k4);
    tt(it)=tt(it-1)+dt;
end;
figure(1)
   plot(tt,x3,'m','linewidth',2)
   xlabel('t');ylabel('x');title('Euler,Feynman,Runge-Kutta')
   pause
   hold off
   
   
   
   
   
   
   
   figure(2)
    plot(tt,abs(x3-xan)./abs(xan),'m',tt,abs(x1-xan)./abs(xan),'r',tt,abs(x2-xan)./abs(xan),'g','linewidth',2)
   xlabel('t');ylabel('x');title('error')
   drawnow