%RSCJ equations to do iv with jj.m
% see jj.m for units
%
global betac Ib
betac= input('What is beta-c?  ');
v=[]; 
I=[];
init = [0 0];
x0 = init';
 for Ib = 0:.1:2;
t0 =0; tf=200;
[t,x] = ode45('jj',t0,tf,x0);
y=x(.1*size(x):size(x),1);
t2=t(.1*size(x):size(x));
v=[v, trapz(t2,y)/(max(t2)-min(t2))];
I=[I, Ib];
x0=[y(size(y))]';
end;
figure;
%v/sqrt(betac) puts V in units of IcRn
plot(v/sqrt(betac),I,'*');
%axis([0 max(I) 0 max(I)]);
title(['Junction with beta_c = ', num2str(betac) ]);
xlabel('voltage in units of IcRn');
ylabel('current in units of Ic');
