%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% damped example: implicit and explicit Euler
% VHS Nov. 2003
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf
global h
global y0
[y1,y2]=meshgrid(-4:0.5:8,-3:0.4:6);
v1=y2;
v2=-3*sin(y1)-0.2*y2;
quiver(y1,y2,v1,v2);
hold on;
k = 0;
h=0.01;
while (k == 0)
    k = waitforbuttonpress;
    point = get(gca,'CurrentPoint');
    point = point(1,1:2);              % extract x and y
    
    % explicit Euler
    z1(1) = point(1);
    z2(1) = point(2);
    for i=1:10/h
        z1(i+1)=z1(i) + h*z2(i);
        z2(i+1)=z2(i) + h*(-3*sin(z1(i))-0.2*z2(i));
    end
    
    % implicit Euler
    y0 = point;
    x1(1) = y0(1);
    x2(1) = y0(2);
    for i=1:10/h
        x0 = y0;
        options = optimset('Display','off');
        [x,fval] = fsolve('imeu',x0,options);
        x1(i+1)=x(1);
        x2(i+1)=x(2);
        y0 = x;
    end
    
    [t,y]=ode23('damped',[0 10],point);
    plot(y(:,1),y(:,2),'r-');
    plot(z1,z2,'g-');
    plot(x1,x2,'b-');
end