function int_flame
%
% Solve IVP:
%   y' = i omega y
%   y(0) = 1
%

y_0 = 1.0E-02;
T = 2/y_0;
n = 100;
h = T/n;
t   = 0;
y = zeros(n+1,1);


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

% exact solution
alpha = 1/y_0 - 1;
y_exact = @(t) 1./(lambertw(alpha*exp(alpha-t))+1);


% plot the solution vs exact solution
figure(1); clf;
fs = 24; %font size

plot(0:h:T,y,'b-',...
     0:h:T,y_exact(0:h:T),'r-',...
     'LineWidth',3,'MarkerSize',10);
xlabel('t','FontSize',fs); ylabel('y(t)','FontSize',fs);

% plot the solution vs exact solution
figure(2); clf;
fs = 24; %font size

plot(0:h:T,y-y_exact(0:h:T)','rx-',...
     'LineWidth',3,'MarkerSize',10);
xlabel('t','FontSize',fs); ylabel('Error','FontSize',fs);

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 y=midstep(t,y,h)
%one step of Midpoint Method
%Input: current time t, current value y,  stepsize h
%Output: approximate solution value at time t+h
k1=ydot(t,y);
k2=ydot(t+h/2,y+h/2*k1);
y =y + h*k2;

function y=trapstep(t,y,h)
%one step of explicit Trapezoid Method
%Input: current time t, current value y,  stepsize h
%Output: approximate solution value at time t+h
k1=ydot(t,y);
k2=ydot(t+h,y+h*k1);
y =y + h/2*(k1+k2);


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 y1=imptrapstep(t,y,h)
% solve the equation y1-h/2*f(t+h,y1)=y+h/2*f(t,y)
% use Newton's method
tol = 1.0E-10;  MaxIt = 100;
y1 = y;
for i=1:MaxIt,
    resid = y1 - h/2*ydot(t+h,y1) - y - h/2*ydot(t,y);
    if abs(resid)/max(abs(y),1) < tol,
        return;
    end
    jac   = 1  - h/2*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
z=y^2-y^3;

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