%E77 Lab2 Task 3
%David Luong, Adem Kader, Mark Piper
%Code adapted from work done by Erik Cheever

% Set Circuit Values
R1=1000;
R2=4000;
R3=1000;
Vt=25E-3;
Is=1E-15;

DT=1E-5;        %Time increment
t=0:DT:.02;    %Time vector

% Allocate memory for arrays.
V1=4*sin(2*pi()*1000*t); %Vin
V2=zeros(size(V1));
V3=zeros(size(V1));


err=1;              %Initialize error to something large to start loop
V0=[0.5 0.5]';     %Initial guess for node voltages [V2 V3]

for i=1:length(V1),
    while err>0.0001,
        e1=exp((V0(1)-V0(2))/Vt);   %Save computation (and typing)
        e2=exp((V0(1)-V1(i))/Vt);   %Save computation (and typing)
        fv1=(V1(i)-V0(1))/R1+(V0(2)-V0(1))/R2-Is*(e2-1)-Is*(e1-1);
        fv2=-V0(2)/R3-(V0(2)-V0(1))/R2+Is*(e1-1);
        dfv1dv1=-1/R1-1/R2-Is*(e2-1)/Vt-Is*(e1-1)/Vt;
        dfv1dv2=1/R2+Is*(e1-1)/Vt;
        dfv2dv1=1/R2+Is*(e1-1)/Vt;
        dfv2dv2=-1/R3-1/R2-Is*(e1-1)/Vt;
        F=[fv1 fv2]';  %Function we wish to find zero of.
        J=[dfv1dv1 dfv1dv2; dfv2dv1 dfv2dv2];  %Jacobian
        Vnew=V0-inv(J)*F;
        err=sum(abs(Vnew-V0)./abs(Vnew+V0)); %error
        V0=Vnew;                             %New guess
    end
    V2(i)=V0(1);           %Store result
    V3(i)=V0(2);
    err=1;                 %Get ready to start loop again.
end

%plot results
plot(t*1000,V1,'r',t*1000,V2,'g',t*1000,V3,'b');
axis([0 2 -4 4]);
grid on
title('Task 3, MatLab');
xlabel('Time (mS)');
ylabel('Voltage (V)');