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

T = 5;
n = 400;
h = T/n;
t = 0;
y = zeros(n+1,1);

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

% plot the trajectory
figure(1); clf;
fs = 24; %font size
xmin = min(real(y));
xmax = max(real(y));
ymin = min(imag(y));
ymax = max(imag(y));

h=plot(real(y),imag(y),'k',...
       real(y(end)),imag(y(end)),'ro',...
       'LineWidth',3,'MarkerSize',10);
xlabel('real(y)','FontSize',fs); ylabel('imag(y)','FontSize',fs);
axis([xmin xmax ymin ymax]); %// freeze axes

animate = 1;

if animate,
    % animate
    for i=1:1:n+1,
        %pause(0.01);
        set(h(1),'XData',real(y(1:i)),'YData',imag(y(1:i)));
        set(h(2),'XData',real(y(i)),  'YData',imag(y(i)));
        drawnow;
    end
end

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 z=ydot(t,y)
%right-hand side of the ODE
omega = 2*pi;
z=j*omega*y;