% stability lecture, June 10th, 2015

close, clear
% Curtis-Hirschfelder example
R = 50;
R = 2000;
f = @(t,y) -R*(y-cos(t));

N = 40; % number of timesteps
%N = 80;
%N = 120;
h = 1.5/N; % step size

t(1) = 0; % initial time
% initial values
w(1) = 0;
y(1) = 0; 
z(1) = 0;


for n=1:N-1
    t(n+1) = t(n) + h;
    % explicit Euler
    %w(n+1) = w(n) + h*f(t(n),w(n));
    % implicit Euler
    y(n+1) = fsolve(@(x) x-y(n)-h*f(t(n+1),x),y(n));
    % trapeziodal rule
    z(n+1) = fsolve(@(x) x-z(n)-h/2*(f(t(n),z(n)) + f(t(n+1),x)),z(n));
end

%hold on
%plot(t,w,'r-d',t,y,'b-o',t,z,'g-*')
plot(t,y,'b-o',t,z,'g-*')
%legend('explicit Euler','implicit Euler', 'TR')
legend('implicit Euler', 'TR')