%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Backward Euler method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for MTH 452-552
%%% this is a simple example how to set up 
%%% a one-step implicit method (here Euler method is used)
%%% to solve a scalar IVP
%%% the data for the problem <f,y0,T,n> can be modified
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% M. Peszynska,  2005/1
%% Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% INSTRUCTIONS: the code can be executed as script or retyped
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% set up the function
f = inline('-5*t.*y.*y+5./t-1/t^2','t','y');
g = inline('h*(-5*t.*y.*y+5./t-1/t^2) + yp - y','t','yp','h','y');
gprime = inline('-h*10*t.*y -1','t','yp','h','y');
phi = inline('1./t','t');
%% fix the initial condition
y0 = 1; 
%% set up the time interval and its discretization
T0 = 1; T  = 5; 
%% stoppping criteria for nonlinear solver
atol = 1e-1; maxit = 20; %% allow at most maxit iterations 
%% the choices below demonstrate three regimes: above, between, below 
%h = .22; choice = 2;  %% BE with Newton method works fine but EE diverges
h = .12; choice = 1;   %% BE with Newton method works fine, EE converges since h is small enough 
%h = .04; choice = 2;  %% very strong constraint - necessary for the fixed pt method to converge
%%
n = floor((T-T0)/h);
%% time steps are stored in t
t = linspace(T0,T,n+1);
%% the numerical solution will be in yi, 
%% for comparison we also show EE solution in ye
ye = zeros(n+1,1) + y0; 
yi = zeros(n+1,1) + y0;
%% initial condition
ye(1) = y0; yi(1) = y0;
maxiter = -inf;
aveiter = 0;
for i=1:n
    %% compute the EE solution
    ye(i+1) = ye(i) + h * f( t(i), ye(i) ); 
    %%% solve by simple iteration the problem  yi(i+1) = yi(i) + h * f( t(i+1), yi(i+1) ); 
    guess = yi(i); tol = inf;
    iter = 0;
    if choice == 1 
        method = 'Picard';
        while tol > atol && iter < maxit
            newguess = yi(i) + h*f(t(i+1),guess);
            %% report for those curious what is the current error
            tol = abs(newguess - yi(i)-h*f(t(i+1),newguess));
            guess = newguess;
            iter = iter + 1;
        end;
    else
        method = 'Newton';
        while tol > atol && iter < maxit
            newguess = guess - g(t(i+1),yi(i),h,guess)/gprime(t(i+1),yi(i),h,guess);
            tol = abs(g(t(i+1),yi(i),h,newguess));
            guess = newguess;
            iter = iter + 1 ;
        end;
    end;
    if iter > maxiter maxiter = iter; end;
    if iter == maxit error (' no convergence for method= ',method); end;
    aveiter = aveiter + iter;
    yi (i+1) = guess;
end;
aveiter = aveiter/n;

fprintf('%s method required iterations: %d max and %g average. Plotting ....',method,...
    maxiter,aveiter);

%% plot the analytical and numerical solution 
plot(t,phi(t),'k',t,ye,'r*');pause,  %%% exact solution is black, EEuler solution is red
plot(t,phi(t),'k',t,yi,'g.');pause,  %%%                          BEuler solution is green
return;


%%% comment out the lines below to exploit addiitonal features of the demo.
%% additional demo: solve the problem using ode45
%[t45,y45]=ode45(f,[1,5],[1]);
%plot(t45,phi(t45),'k',t45,y45,'b');
return;
%% plot the analytical solution and then plot individual steps of numerical
%% solution
plot(t,phi(t),'k'); hold on;
for i=1:n 
    plot([t(i),t(i+1)],[ye(i),ye(i+1)],'r*',[t(i),t(i+1)],[yi(i),yi(i+1)],'g.');pause,
end;
hold off;
