clear all;
global a;%全局變量的設置至關重要，是整個程序的關鍵
a=1;
BB=[];
for a=[1:0.5:9]
    % 選用二三階的龍格庫塔函數，其中0到5的步進為x的取值范圍，pi/6為x取值為0時y的值
    [x,y]=ode23(@fangcheng,[0:0.001:5],[pi/6]);
    b=find(abs(y-pi/3)<1e-3);%尋找y=pi/3的元素標號
    I1=a*5-5*sin(x(b)/4.8*0.144*pi);%根據相關公式計算電流值
    [x,y]=ode23(@fangcheng,[10:0.001:15],[pi/6]);%二次取值
    c=find(abs(y-pi/3)<1e-3);
    I2=a*5-5*sin(x(c)/4.8*0.144*pi);
    s=length(I1);s1=length(I2);
    A=zeros(1,s+s1);
    A(1:s)=I1';
    A((s+1):(s+s1))=I2';
    B=zeros(2,s+s1);%創建數組
    B(1,:)=A;
    V=(a-0.866)*10;
    B(2,:)=V';
    for r=[1:s+s1]%篩選出去一部分數值
        if B(1,r)<5
            B(2,r)=0;
        end
    end
    BC=[];
    for i=1:length(B)%將數組中的列重排
        BC(:,i)=B(:,(length(B)-i+1));
    end
    BB=[BB BC];%將所有結果存入大數組
end
figure(1);
stairs(BB(2,:),BB(1,:));%階梯畫圖
title('微波輻照下約瑟夫森結的I-V特性');
xlabel('\fontsize{10}\bfV(mV)\rightarrow')
ylabel('\fontsize{10}\bfI(mA)\rightarrow')