%
clear all;
%
blank=sprintf(' ');
disp(blank);
disp('   Duffing_pendulum.m');
disp(blank);
disp('   by Tom Irvine');
disp('   Email: tomirvine@aol.com');
disp(blank);
%
out1=sprintf(' Enter omegan (rad/sec) ');
omegan=input(out1);
%
out1=sprintf(' Enter omega (rad/sec) ');
omega=input(out1);
%
a=(1/8)*omegan^2;
c=-omegan^2+omega^2;
%
ia=1;
ib=1;
ic=1;
%
for(i=1:10000)
    df=0.001;
    F=(i-1)*df;
    if(F >= 0.35 ), break, end    
%
    clear p;
    clear r;
    p=[a 0 c F]; 
    r=roots(p);  
%
    iflag=0;
%
    if( imag(r(1))==0 & imag(r(2))==0 & imag(r(3))==0)
        dd(1)=abs(r(1));
        dd(2)=abs(r(2));   
        dd(3)=abs(r(3));
        ds=sort(dd);
        disp1(i)=ds(1);
        disp2(i)=ds(2);
        disp3(i)=ds(3);  
        force1(ia)=F;
        force2(ib)=F;
        force3(ic)=F;
        ia=ia+1;
        ib=ib+1;
        ic=ic+1;
        if(disp1(i) >= pi ), break, end
        if(disp2(i) >= pi ), break, end  
        if(disp3(i) >= pi ), break, end
        iflag=1;
    end
%
   if( imag(r(1))==0 & imag(r(2))~=0 & imag(r(3))~=0)
       disp1(i)=abs(r(1)); 
       force1(ia)=F;    
       ia=ia+1;  
   end  
   if( imag(r(1))~=0 & imag(r(2))~=0 & imag(r(3))==0)
       disp1(i)=abs(r(3));    
       force1(ia)=F;
       ia=ia+1;  
   end   
%    
%    out1=sprintf(' %f %f %f',F,disp1(i),disp2(i));
%    disp(out1);
end
%
if( iflag==1 )
   plot(force1,disp1,force2,disp2,'-.',force3,disp3,'--');
else
   plot(force1,disp1)
end
%
out5 = sprintf(' omegan=%g rad/sec   omega=%g rad/sec',omegan,omega);
title(out5);
%
xlabel(' Scaled Moment [M/(m L^2) ]  (rad/sec^2)');
ylabel(' Angular Displacement (rad) ');
grid on;