% solves y'=F(x=t,y) and loops over size of time step
% plots error as function of time step

clear 
format long
y0=1;
  
dt0=0.1;
tmax=.995;
  
jmax=14;
dtfactor=2;

for j=1:jmax
  dt=dt0/(dtfactor^(j-1))
  dtj(j)=dt;
end

for j=1:jmax
 dt=dtj(j);
 numsteps=floor(tmax/dt)+1;
 x(1:numsteps)=0;
 y(1:size(y0),1:numsteps)=0; 

 y(1,1)=y0;

 cput1=cputime;

 for k=1:numsteps

  if x(k)+dt > tmax
   dt=tmax-x(k);
  end

% Forward Euler
  [Fy] = F(x(k),y(:,k));
  yk1 = dt*Fy';
  y(:,k+1)=y(:,k)+yk1;
  x(k+1)=x(k)+dt;

 end

 elapsed_time=cputime-cput1

 yexact=sqrt(1-x);
 solution=y(:,numsteps+1)
 solex=yexact(numsteps+1)
 error(j)=abs(y(1,numsteps)-yexact(numsteps));

 figure(1);
  plot(x(:),y(1,:),'b+');
  xlabel('Time')
  ylabel('y')
  hold on
  plot(x(:),yexact(:),'r*');
  hold off

 figure(2)
  loglog(dtj(1:j),error(1:j),'-ro','LineWidth',2,...
         'MarkerSize',14,'MarkerFaceColor',[1,0,0]);
  xlabel('Time Step dt')
  ylabel('Error')
end
