% Code from A Brief History of MATLAB - MATLAB & Simulink:
%

clear
mu = 5;
vdp = @(t,y) [y(2); mu*(1-y(1)^2)*y(2)-y(1)];
tspan = [0 30];
y0 = [0 0.01]';
[t,y] = ode23s(vdp,tspan,y0);
plot(t,y)
legend({'y','dy/dt'})
xlabel('t')
ylabel('y')
title('Moler''s Van der Pol Oscillator')

%

clear
omega = 1;
shm = @(t,y) [y(2); -omega^2*y(1)];
tspan = [0 30];
y0 = [0 0.01]';
[t,y] = ode23s(shm,tspan,y0);
plot(t,y)
legend({'y','dy/dt'})
xlabel('t')
ylabel('y')
title('Simple Harmonic Motion')

figure
plot(y(:,1),y(:,2))
xlabel('x')
ylabel('y')
title('SHM Trajectory')
axis equal

%

y0 = 0:.1:2;
t0 = 0;
tfinal = 15;
for k = 1:length(y0)
    [t,p] = ode45(shm,[t0 tfinal],[0 y0(k)]);
    plot(p(:,1),p(:,2),'b')
    hold on
end
title('Phase Plot for SHM')
xlabel('x')
ylabel('t')
hold off
axis equal

%

clear
omega = 1;
nlp = @(t,y) [y(2); -omega^2*sin(y(1))];
tspan = [0 30];
y0 = [0 0.01]';
[t,y] = ode23s(nlp,tspan,y0);
plot(t,y)
legend({'y','dy/dt'})
xlabel('t')
ylabel('y')
title('Nonlinear Pendulum')

x0 = -6:1:6;
y0 = -pi:.5:pi;
t0 = 0;
tfinal = 5;
for j = 1:length(x0)
    for k = 1:length(y0)
        [t,p] = ode45(nlp,[t0 tfinal],[x0(j) y0(k)]);
        plot(p(:,1),p(:,2),'b')
    end
    hold on
end
title('Phase Plot for Nonlinear Pendulum')
xlabel('x')
ylabel('y')
hold off
axis([-3*pi,3*pi,-pi,pi])