clear
T=20;%final time
tspan=[0 T];%time interval
y0=[2;0]; %initial condition
opts = odeset('RelTol',1e-12,'AbsTol',1e-12);
%vdp1.m is MATLAB built-in van der Pol ODEs for mu = 1
[t,y] = ode45(@vdp1,tspan,y0,opts);%variable step size non-stiff ODE solver
%[t,y] = ode15s(@vdp1,tspan,y0);%variable step size stiff ODE solver

%plot the solutions, notice the time points are not uniform
plot(t,y(:,1),'-o',t,y(:,2),'--x')
title('Solution of van der Pol Equation ($\mu = 1$) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1(t)','y_2(t)')
set(gcf, 'Position', get(0, 'Screensize')/2);
legend('-DynamicLegend','Location', 'Best');

%solve Lorenz equation
tspan=[0 100];
y0=[1;1;1];
[t,y]=ode15s(@Lorenz,tspan,y0);
%Figure 6.17
plot(y(:,1),y(:,3))

% [t,y] = ode15s(@vdp1000,[0 3000],[2 0]);
% plot(t,y(:,1),'-o');
% hold on
% [t,y2] = ode45(@vdp1000,[0 3000],[2 0]);
% plot(t,y2(:,1),'--s');