%  This is a Matlab script that uses the Euler Method 
%  (and also the Improved 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.  There is also an exact solution calculated
%  and plotted but the appropriate exact solution depends on the 
%  right hand side defined in feuler.
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
%
%
% 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');