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

global omega

omega = 2*pi;

y_exact = @(t) exp(j*omega*t);

T = 5;
n = 200;

% arrays for storing h and error
H = [];
E = [];

% perform mesh refinement
for k=1:12,
    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
    
    % compute the error
    g = abs(y-y_exact(h*(0:n)'));
    E = [E,max(g)];
    H = [H,h];
    
    % refine the mesh
    n = n*2;
end

% plot the max global error vs h
figure(1); clf;
fs = 24; %font size
loglog(H,E,'kx-',...
       H,H.^2,'r-',...
       'LineWidth',3,'MarkerSize',10);
xlabel('h','FontSize',fs); ylabel('error','FontSize',fs);
grid on;
h=legend('Global error','Asymptotics');
set(h,'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 z=ydot(t,y)
%right-hand side of the ODE
global omega

z=j*omega*y;