% Implicit vs explicit Euler method for stiff problems
function oppgave_6_6_2

y_0 = 0.5;
a = 0;
b = 20;
n = 66;
h = (b-a)/n;
t = a;
ye = zeros(n+1,1);
yi = zeros(n+1,1);

% start the integration loop
ye(1) = y_0; yi(1) = y_0;
for i=1:n,
    ye(i+1) = eulerstep(t,ye(i),h);
    yi(i+1) = impeulerstep(t,yi(i),h);
    t = t+h;
end

% plot the explicit vs implicit approximation
figure(1); clf;
fs = 24; %font size

plot(a:h:b,yi,'b-',...
     a:h:b,ye,'r-',...
     'LineWidth',3,'MarkerSize',10);
xlabel('t','FontSize',fs); ylabel('y(t)','FontSize',fs);
title(sprintf('Explicit vs implicit Euler, h=%e',h));
legend('Explicit Euler','Implicit Euler');


function y=eulerstep(t,y,h)
%one step of forward Euler Method
%Input: current time t, current value y,  stepsize h
%Output: approximate solution value at time t+h
k1=ydot(t,y);
y =y + h*k1;

function y1=impeulerstep(t,y,h)
% solve the equation y1-h*f(t+h,y1)=y
% use Newton's method
tol = 1.0E-10;  MaxIt = 100;
y1 = y;
for i=1:MaxIt,
    resid = y1 - h*ydot(t+h,y1) - y;
    if abs(resid)/max(abs(y),1) < tol,
        return;
    end
    jac   = 1  - h*ydot_prime(t+h,y1);
    y1    = y1 - resid/jac;
end
fprintf('Warning: Newton''s iteration did not converge at t=%e\n',t);


function z=ydot(t,y)
%right-hand side of the ODE
% oppgave 6.6.2 (a)
z=6*y-3*y^2;
% oppgave 6.6.2 (b)
%z = 10*y^3-10*y^4;

function z=ydot_prime(t,y)
%derivative of the right-hand side with respect to y
% (used for implicit methods)
% oppgave 6.6.2 (a)
z=6-6*y;
% oppgave 6.6.2 (b)
%z = 30*y^2-40*y^3;