%                        EulerSample.m
%
% A sample Matlab script that implements the Euler's method for
% the ODE
%
% dy/dt = F(y,t)
%
% where the function F(y,t) is specified by the m-file F.m
%
%method=1 for forward Euler, 2 for backward Euler, 3 for trapezoidal method
%
tInitial    = 0.0;                      % Initial time
tFinal      = pi;                      % Final time
yInitial    = 0.0;                      % Initial value of y
N           = 150
dt          = (tFinal - tInitial)/N;    % Timestep determination
y           = zeros(N+1,1);             % Arrays to hold solution values
t           = zeros(N+1,1);
error           = zeros(N+1,1);
yex           = zeros(N+1,1);             % Arrays to hold solution values
lambda      =-100;
title1='Stiff ODE Example';
%
%////////////////////////////////////////////////////////
%     Computing the solution using Euler's method
%////////////////////////////////////////////////////////
%

for method=1:2
    t(1)  = tInitial;
    y(1)  = yInitial;
    yex(1)=yInitial;
    error(1)=0;

    if method==1
        for(i = 1:N)
            t(i+1) = t(i) + dt;
            y(i+1) = y(i) + dt*lambda*(y(i)-G(t(i),lambda));
        end
    elseif method==2
        for(i = 1:N)
            t(i+1) = t(i) + dt;
            y(i+1) = BEInv(y(i),t(i),lambda,dt);
        end
    elseif method==3
        for(i = 1:N)
            t(i+1) = t(i) + dt;
            %         y(i+1) = y(i) + dt*F(y(i),t(i));
            %         y(i+1) = y(i) + dt*FHW1(y(i),t(i));
            y(i+1) = ((2-100*dt)*y(i) + dt*(100*sin(t(i+1))+100*sin(t(i))+cos(t(i+1))+cos(t(i))))/(2+100*dt);
            %error(i+1)=y(i+1)-sin(t(i+1));
        end
    end
    for(i = 1:N)
        yex(i+1)=sin(t(i+1))+yInitial*exp(lambda*t(i+1));
        error(i+1)=y(i+1)-yex(i+1);
    end

    % Now assign names to y,error
    if method==1
        yF=y;errorF=error;
    elseif method==2
        yB=y;errorB=error;
    elseif method==3
        yT=y;errorT=error;
    end
end %end of loop over method
%
%    Set plot limits
%
tPlotMin =     tInitial;
tPlotMax =     tFinal;
yPlotMin =     0.0;
yPlotMax =     2.0;

% plot,
figure(1)
plot(t,yF,'-r',t,yB,'-b',t,yex,'-kx');
xlabel('t');ylabel('y');
title([title1 ' dt= ' num2str(dt) ', N= ' num2str(N)])
legend('Forward Euler','Backward Euler','Exact')

figure(2)
semilogy(t,abs(errorF),'-r',t,abs(errorB),'-b');
xlabel('t');ylabel('Error');
%title(['Forward Euler. Stiff example dt= ' num2str(dt)])
title([title1 ' dt= ' num2str(dt) ', N= ' num2str(N)])
legend('Forward Euler','Backward Euler','Exact')