clear
% solve y' = y*sin(t) using Euler's method
% setting up grid in time
for j=1:5
    n=10*2^j;
t=linspace(0,1,n);
for i=1:n
  t(i) = 2*(i-1)/(n-1);
end
dt=t(2)-t(1);
% 
y(1) = 1.0
% Euler's method
% for i=1:n-1
%   y(i+1) = y(i) + dt*(y(i)*sin(t(i)));
% end
% Improved Euler's method
%  for i=1:n-1
%      y1=y(i);
%      t1=t(i);
%      f1= y1*sin(t1);
%      y2=y1+dt*f1;
%      t2=t1+dt;
%      f2= y2*(sin(t2));
%     y(i+1) = y(i) + 0.5*dt*(f1 + f2);
% end
%RK4
 for i=1:n-1
    k1= dt*y(i)*sin(t(i));
    k2= dt*(y(i)+k1/2)*sin(t(i)+dt/2);
    k3= dt*(y(i)+k2/2)*sin(t(i)+dt/2);
    k4= dt*(y(i)+k3)*sin(t(i)+dt);
    y(i+1) = y(i) + (k1+2*k2+2*k3+k4)/6.0;
end

   
% error analysis
error(1)=0.0;
for i=1:n-1
   c = exp(1);
   exact(i+1) = c*exp(-cos(t(i)+dt));
   error(i+1)=abs(y(i+1)-exact(i+1));
end
norm(error,'inf')
%plot(t,y,'r')
%plot(t,error,'g')
hold on
end
plot(t,exact,'r')


