% solves ode with Euler or RK4 with
% time-step control and Richardson extrapolation

clear
format long
cputime1=cputime;

global diff;

dtfactor=2;

% for tmax=0.995 Euler without extrapol goes for jmax=13 to error=3e-4
% (with 300,000 steps)
% for tmax=0.995 Euler with extrapol goes for jmax=13 to error=2e-8
% (with 300,000 steps)
% for tmax=0.995 RK4  without extrapol goes for jmax=9 to error=1e-10
% (with 800 steps)
% for tmax=0.995 RK4  with extrapol goes for jmax=9 to error=1e-12
% (with 800 steps)

% Forward Euler
jmax=13;
tol0=.1;
p=2;

% RK4
jmax=20;
tol0=.1;
p=5;

% allocate memory to the arrays to speed up the code
ntmax=1000000;
% sol stores position of particle for all times
sol=zeros(ntmax,1);
% for variable time step need to keep track of deltat and t separately
deltat=zeros(ntmax,1);
t=zeros(ntmax,1);

tmax=.999;

for j=1:jmax
  tolj(j)=tol0/(dtfactor^((2-1)*(j-1)))
end

for j=1:jmax
tol=tolj(j)
y=[1];
time=0;
dt=0.01;
nt=1;
sol(1)=y(1);
t(1)=time;
deltat(1)=dt;

while time < tmax
 
 timep=time+dt;

 if timep > tmax
   dt;
   time;
   dt=tmax-time;
 end

%  [time,y,dt,success]=steprk(time,y,dt,tol);
  [time,y,dt,success]=steprkopt(time,y,dt,tol,tmax,p);
%  [time,y,dt,success]=steprkoc(time,y,dt,tol);
%  [time,y,dt,success]=steprkocwe(time,y,dt,tol);
%  [time,y,dt,success]=stepeuleroc(time,y,dt,tol);
%  [time,y,dt,success]=stepeuler(time,y,dt,tol);
%  [time,y,dt,success]=stepeuleropt(time,y,dt,tol,tmax,p);
%   [time,y,dt,success]=stepeuleropt(time,y,dt,tol,tmax);
  if success > 0
   nt=nt+1;
   sol(nt)=y(1);
   t(nt)=time;
   deltat(nt)=dt;
  end
  if mod(nt,1000)==0 
    nt 
  end

end
 
ntj(j)=nt;
nt
time
 
ufinal=sol(nt)
uexact=sqrt(sqrt(1-time));
numerror(j)=abs(ufinal-uexact);
deltat2=deltat(floor(nt/4):nt-1);
dtmax(j)=max(deltat2);
dtmin(j)=min(deltat2);

cputime2=cputime;
cputimetotal=cputime2-cputime1;
cputime1=cputime2;

%figure(1);
%plot(t(1:nt),sol(1:nt));
%ylabel('u');
%figure(2);
%plot(t(1:nt),deltat(1:nt),'+');
%ylabel('Delta t');
%xlabel('t')
figure(3);
semilogy(t(1:nt),deltat(1:nt),'+');
ylabel('Delta t');
xlabel('t')
figure(4)
loglog(dtmax(1:j),numerror(1:j),'-bo','LineWidth',2,...
         'MarkerSize',14,'MarkerFaceColor',[0,0,1])
hold on
loglog(dtmin(1:j),numerror(1:j),'-rs','LineWidth',2,...
         'MarkerSize',14,'MarkerFaceColor',[1,0,0])
ylabel('error')
xlabel('dt')
hold off
figure(5)
loglog(ntj(1:j),numerror(1:j),'-ro','LineWidth',2,...
         'MarkerSize',14,'MarkerFaceColor',[0,0,1])
xlabel('nt')
ylabel('error')
figure(6)
semilogx(tolj(1:j),numerror(1:j)./tolj(1:j),'ro','LineWidth',2,...
         'MarkerSize',14,'MarkerFaceColor',[0,0,1])
	 xlabel('tolerance')
	 ylabel('error/tolerance')
end

