%Math 316 - Spring 2017
%Numerical ODEs
%MATLAB code for EXAMPLE 1 
%Written by M. Motamed @ UNM 
%===============================================================
%Part I - plot the exact solution and three numerical solutions 
%===============================================================
clear all

t=linspace(0,1,1000);
yexact=3*exp(0.5*t.^2)-t.^2-2;

figure
plot(t,yexact,'b-')
xlabel('t')
ylabel('y(t)')
set(gca,'fontsize', 22);
hold on 


N=(1:3);
H=0.2./(2.^N); % to generate a decreasing sequence of step sizes 
color={'ro-' 'k*-' 'gv-'};

for k=1:length(H)
   h=H(k);           %step size
   t=(0:h:1);
   y=zeros(size(t)); %to pre-allocate a vector for the solution
   y(1)=1;           %to set the initial solution at t=0
   for i=1:length(t)-1
     f=t(i)*y(i)+t(i)^3;
     y(i+1)=y(i)+h*f;
   end
   plot(t,y,color{k})
end

legend('exact solution',['Euler (h=' num2str(H(1)) ')'],'Euler (h=1/10)','Euler (h=1/20)')

%===============================================================
%Part II - plot the numerical error at t=1
%===============================================================

N=(1:15);
H=0.2./(2.^N);        % to generate a decreasing sequence of step sizes 
E=zeros(size(H));     %to pre-allocate a vector for the error
yexact1=yexact(end);  %exact solution at t=1

for k=1:length(H)
  h=H(k);         %step size
  n=round(1/h);   %number of sub-intervals
  t=0;
  y=1; 
  for i=1:n
    f=t*y+t^3;
    y=y+h*f;
    t=t+h;
  end
  E(k)=abs(yexact1-y);
end

figure
loglog(H,E,'o-')
hold on
loglog(H,H,'r-')
xlabel('h')
ylabel('error')
set(gca,'fontsize', 22); 
legend('Euler error at t==1','h')

