%  This is a Matlab script that uses the Backward 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=100;
DELTAT=(TF-T0)/NK;
T=T0;Y=Y0;
tk=T0;yk=Y0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Backward Euler Method (nonlinear solve via secant method)
%
%%%%%%%
%
% define tolerance, max iteration and initial guess for secant method
%
sectol=DELTAT*10^(-2);
itmax=20;
y1=Y0;
y2=Y0+DELTAT;         % a somewhat uneducated guess for y2 ...
%
% execute the integration loop
%
for k=1:NK;
tk=tk+DELTAT;
yold=yk;
f1=y1-yold-DELTAT*feuler(tk,y1);
f2=y2-yold-DELTAT*feuler(tk,y2);
%
% set or reset parameters for secant method solve
%
dy=1;
i=1;
%
% implement a Secant Method to solve y(i+1)-y(i)-h*f(t(i+1),y(i+1))=0;
%
    while (dy>sectol & i<itmax)
        fp=(f2-f1)/(y2-y1);
        yk=yk-f2/fp;
        y1=y2;
        f1=f2;
        y2=yk;
        f2=y2-yold-DELTAT*feuler(tk,y2);
        dy=abs(-f2/fp);
        i=i+1;
    end
%
    if (i==itmax)
       disp('itmax exceeding in nonlinear solve'),i
    end
%
    T=[T tk];
    Y=[Y 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)      %normBE
%
% Plot the solutions
%
subplot(2,1,1);        % upper 
   hold on;
   plot(T,Y,'o-');     % Euler Solution
   plot(T,YEXACT);
   axis([min(T) max(T) min(YEXACT) max(YEXACT)])
   grid;
   xlabel('t-axis');
   ylabel('y-axis');
   title('Backward Euler (o), 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');