%RSCJ equations
%
global betac Ib Rsg
betac= input('What is beta-c?  ');
Rsg=input('Ratio of Rsg/Rn  ?');
Ib = input('What is the current?  ');
init = input('initial starting point   ');
t0 =0; tf=200;
x0 = init';
[t,x] = ode45('jjnlinb',t0,tf,x0);
figure;
plot(t,x(:,1));
title(['voltage vs time']);
pause;
figure;
plot(t,sin(x(:,2)));
title(['phase vs time']);
pause;
%figure;
%plot(t,(x(:,1)));
%
%to have the phase portrait
%pause;
figure;
plot(rem(x(:,2),2*pi),x(:,1),'.');
%axis([-4 4 -1 1]);
title(['Junction with beta_c = ', num2str(betac), 'and I= ', num2str(Ib)]);
xlabel('theta');
ylabel('theta-dot');
