clc;
clear;
% 此題為龍格-庫塔法算法的應用


A      = 10;
B      = 1;

para.A = A;
para.B = B;
tic;
[y t] = Gill(@vdp_3,[1 9],[1 5],[.01,1e-2,10],para);
Time  = toc;
figure
plot(t,y(1,:),'r-',t,y(2,:),'b:');
xlabel('Time s ');

legend('y1','y2');
title(['耗時:' num2str(Time)]);