%  This is a Matlab script that uses the Euler Method to 
%  solve the first order differential equation y'=f(t,y),
%  and plots the solution from TO to TF.  The functions 
%  to be evaluated are given in a separate function m-file,
%  called feuler(t,y).  The first part of this file solves
%  the equation numerically.  The second part plots the 
%  direction fields.
hold off
clear all
T0=0;
Y0=1;
TF=1;
NK=200;
DELTAT=(TF-T0)/NK;
T=T0;Y=Y0;
tk=T0;yk=Y0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Euler Method
%
for k=1:NK;
    k1=feuler(tk,yk);
    yk=yk+DELTAT*k1;
    tk=tk+DELTAT;
    T=[T tk];
    Y=[Y yk];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Improved Euler Method (aka Huen's Method)
%
TI=T0;YI=Y0;
tk=T0;yk=Y0;
for k=1:NK;
    k1=feuler(tk,yk);
    k2=feuler(tk+DELTAT,yk+DELTAT*k1);
    yk=yk+0.5*(k1+k2)*DELTAT;
    tk=tk+DELTAT;
    TI=[TI tk];
    YI=[YI yk];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Exact Solution
%
%%%YEXACT=exp(T).*sin(5*T);
%%%YEXACT=T./(1+T.^2);
YEXACT=1./(1+90*T).^(1/3);
norm(Y-YEXACT)/norm(YEXACT)      %normE
norm(YI-YEXACT)/norm(YEXACT)     %normIE
%
% results show order one accuracy of Euler
%  [reducing DELTAT by factor of 2 reduces error by factor of 2]
% results show second order accuracy of Improved Euler
%  [reducing DELTAT by factor of 2 reduces error by factor of 4]
%    [for y'=y+5cos(5t)e^t : YEXACT = (e^t)*sin(5t)]
%
%NK   DELTAT              normE           normIE
%5    0.2            0.75009326741153    0.08975239661325
%10   0.1            0.42752738986351    0.02269916540100
%20   0.05           0.22975268766778    0.00572730133133
%40   0.025          0.11934143704239    0.00143948140520
%80   0.0125         0.06085405033235    3.608907620717458e-04
%160  0.00625        0.03073185301058    9.035403879238091e-05
%320  0.003125       0.01544329893073    2.260516006677027e-05
%640  0.0015625      0.00774114091936    5.653385637528957e-06
%1280 7.8125e-04     0.00387546203620    1.413609176290171e-06
%
% results show order one accuracy of Euler
% and second order accuracy of Improved Euler
%     [for y'=1/(1+t^2) - 2y^2 : YEXACT=t/(1+t^2)]
%
%NK    DELTAT                 normE               normIE
%25    0.04          1.978642483772563e-002   7.198195290059719e-004
%50    0.02          9.857599507246146e-003   1.763005471500330e-004  
%100   0.01          4.920149008542324e-003   4.362833063190092e-005  
%200   0.005         2.457942615868702e-003   1.085183503841653e-005   
%400   0.0025        1.228442138265047e-003   2.706089848448282e-006  
%800   0.00125       6.140892490769361e-004   6.756661375833888e-007  
%1600  0.000625      3.070117283889323e-004   1.688096296240505e-007  
%3200  0.0003125     1.534976474940130e-004   4.218905905586938e-008  
%
% results show order one accuracy of Euler
% and second order accuracy of Improved Euler
%     [for y'=-30*y^4 : YEXACT = (1+90t)^(-1/3)]
%
%NK    DELTAT                 normE               normIE
%25    0.04            1.446475234464752e+00   1.491060321350094e-01
%50    0.02            2.016793104847440e-01   7.370273641098403e-03
%100   0.01            5.673729019638415e-02   3.875167766364247e-03
%200   0.005           2.288144770131952e-02   1.355518370487168e-03
%400   0.0025          1.061958979902291e-02   3.489992448877024e-04
%800   0.00125         5.153790709417700e-03   8.632757679926444e-05
%1600  0.000625        2.542534012885682e-03   2.136434319514973e-05
%3200  0.0003125       1.263170609717669e-03   5.308789742970363e-06
%
%
% Plot the solutions
%
subplot(2,1,1);        % upper 
   hold on;
   plot(T,Y,'o-');     % Euler Solution
   plot(TI,YI,'x-');   % Improved Euler
   plot(T,YEXACT);
   axis([min(T) max(T) min(YEXACT) max(YEXACT)])
   grid;
   xlabel('t-axis');
   ylabel('y-axis');
   title('Euler (o), Improved Euler (x), Exact');
%
% This is the part of the file that computes direction
% fields associated with the equation y'=f(t,y).
%
TMIN=T0;
TMAX=TF;
DTPLOT=0.1*(TF-T0);
YMIN=min(YEXACT);
YMAX=max(YEXACT);
DYPLOT=0.1*(YMAX-YMIN);
tj=[TMIN:DTPLOT:TMAX];     % define the time domain for plotting
yj=[YMIN:DYPLOT:YMAX];     % define the y domain for plotting
[TM,YM]=meshgrid(tj,yj);   % creates matrices T and Y containing grid 
                           % coordinates as entries
%
% define rhs functions of your choice (called F1) below
% note that F1 is a matrix
%
%%%F1=YM+5.*cos(5.*TM).*exp(TM);
%%%F1=1./(1+TM.^2)-2*YM.^2;
F1=-30*YM.^4;
%
%
[nn,mm]=size(F1);        % this line identifies the size of the matrix F1
M1=ones(nn,mm);          % this line defines the matrix M1 to be a matrix
                         % of all ones and of the same size as matrix F1
M1=M1./sqrt(1+F1.^2);    % provide normalization of arrows to unit length
F1=F1./sqrt(1+F1.^2);    % provide normalization of arrows to unit length
s=0.30;                  % choice of scale factor for quiver command
                         % choose big or small to adjust length of arrows
subplot(2,1,2);   % lower
   hold on;
   axis([TMIN TMAX YMIN YMAX]);
   quiver(TM,YM,M1,F1,s);
   plot(T,YEXACT);
   grid;
   xlabel('time'); ylabel('y');title('exact solution and direction fields');
