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