Trying to use Newton Raphson Method to solve circuit design

I need help to diagnose why my code doesn't seem to be working probably. I have a feeling it mught be something mundane, but I keep banging my head with no results. Here is the assignment.:

In Task B the circuit in Figure 2 will be implemented in MATLAB. Again, in this case we will assume that the input signal is represented by the AC excitation given in Equation (11). 𝑣 𝑡 can be approximated by following the steps:

1. Sample the time variable over one cycle of time as: W = 1; T = 2*pi/W; tn = 0:T/40:T;

2. For each time sample, use the Newton-Raphson method to find the solutions of 𝑣 , and 𝑣 , . I.e., Solve the system:

3. The output samples (𝑣 , ) can now be calculated using either Equation (6.a) or Equation (6.b). Note that the vectors 𝑣 , and 𝑣 , must have the same size.

4. Once you implement this algorithm generate the following plots:

a. Plot 𝑣 vs. 𝑣 . Use proper labels for the axes in your graph.

b. Create a plot that includes two curves: the first one corresponds to 𝑣 vs. t, and the second one shows 𝑣 vs. t. Use proper labels for the axes in your graph and add a legend.

Here is the code for the Newton Raphson algorithim

clear
clc
w = 1;
T = 2*pi/w;
tn = 0:T/40:T;
I = 10^-15;
Vt = 25.9*10^-3;
Vo = 10;
V1 = 6;
V2 = 6;
R = 50*10^3;

x(1) = 0.7; %Initial Value X
y(1) = 0.7; %Initial Value Y

e_d = 1*10^(-9)*100; %Desired Error

f1 = @(x,y) I.*(exp(x(i)./Vt)-exp(y(i)./Vt)-(Vo.*cos(tn).*V1.*x(i))./R);
f1dx = @(x) (I/Vt).*(exp(x(i)./Vt)-(1./R));
f1dy = @(y) -(I/Vt).*(exp(y(i)./Vt)-(1./R));
f2 = @(x,y) -Vo.*cos(tn).*I.*(exp(x(i)./Vt)-exp(y(i)./Vt))+(R.*I.^2)*(exp(x(i)./Vt)-exp(y(i)./Vt).^2)+(x(i)+V1).*I.*(exp(x(i)./Vt)-1)+(y(i)+V2).*I.*(exp(y(i)./Vt)-1);
f2dx = @(x)  ((-Vo.*cos(tn).*I)./Vt)*exp(x(i)./Vt)+(2.*R.*I.^2)*(exp(x(i)./Vt)-exp(y(i)./Vt).*(exp(x(i)./Vt))+(I.*(((x(i)+V1)./Vt)*(exp(x(i)./Vt))+(exp(x(i)./Vt)-1));
f2dy = @(y) ((Vo.*cos(tn).*I)./Vt)*exp(y(i)./Vt)+(2.*R.*I.^2)*(exp(x(i)./Vt)-exp(y(i)./Vt).*(exp(y(i)./Vt))+(I.*(((y(i)+V2)./Vt)*(exp(y(i)./Vt))+(exp(y(i)./Vt)-1));

J = [f1dx f1dy;f2dx f2dy];

i = 1; %Loop Index - x
j = 1; %Loop Index - y
while (1)
    [x(i+1),y(i+1)] = [x(i);y(i)]-inv(J).*[f1;f2]; %Netwon-Raphson general formula
    %x(i+1)=x(i) - f1(x(i))./f1dx(x(i)); %Newton-Raphson Formula 1
    %y(i+1)=y(i) - f2(y(i))./f2d(y(i)); %Newton-Raphson Formula 2
    if i>1 && j>1
        e_A(i-1) = abs((x(i+1)-x(i))./x(i+1))*100; %Approx.  error x(%)
        e_A(j-1) = abs((y(i+1)-y(i))./y(i+1))*100; %Approx. error y(%)
        if (abs(e_A(i-1)<e_d) && abs(e_A(j-1)<e_d))
            break %If approx. error < desired error terminate loop
        end
    end
    i=i+1; %Update count for next iteration
    Vin= Vo*cos(tn);
    Vout = V1.*eye(i)+x(i); 
end

%Plot Vout vs. Vin 
plot(Vout,Vin)
 xlabel("Vout"), ylabel("Vin"), legend("\Vout","\Vin"), title("Vout v Vin")

%Multi-Plot
fplot(tn, Vout)
 xlabel("t"),ylabel("Vout"), legend("\t","\Vout"), title("Vout v t")
hold on
plot(tn, Vin)
 xlabel("t"),ylabel("Vin"), legend("\t","\Vin"), title("Vin v t")

Answer 1

[Torsten]

w = 1;
T = 2*pi/w;
tn = 0:T/40:T;
I = 1e-15;
Vt = 25.9e-3;
Vo = 10;
V1 = 6;
V2 = 6;
R = 50e3;
x0 = [0.7;0.7]; 
e_d = 1e-9; 
for i = 1:numel(tn)
   t = tn(i);
   f1 = @(x) I*(exp(x(1)/Vt)-exp(x(2)/Vt))-(Vo*cos(t)-V1-x(1))/R;
   f2 = @(x) -Vo*cos(t)*I*(exp(x(1)/Vt)-exp(x(2)/Vt))+R*I^2*(exp(x(1)/Vt)-exp(x(2)/Vt))^2+...
             (x(1)+V1)*I*(exp(x(1)/Vt)-1)+(x(2)+V2)*I*(exp(x(2)/Vt)-1);
   f = @(x)[f1(x);f2(x)];
   df1dx1 = @(x)I*exp(x(1)/Vt)/Vt+1/R;
   df1dx2 = @(x)-I*exp(x(2)/Vt)/Vt;
   df2dx1 = @(x)-Vo*cos(t)*I*exp(x(1)/Vt)/Vt+2*R*I^2*(exp(x(1)/Vt)-exp(x(2)/Vt))*exp(x(1)/Vt)/Vt+...
               I*(exp(x(1)/Vt)-1)+(x(1)+V1)*I*exp(x(1)/Vt)/Vt;
   df2dx2 = @(x)Vo*cos(t)*I*exp(x(2)/Vt)/Vt+2*R*I^2*(exp(x(1)/Vt)-exp(x(2)/Vt))*(-exp(x(2)/Vt))/Vt+...
               I*(exp(x(2)/Vt)-1)+(x(2)+V2)*I*exp(x(2)/Vt)/Vt;
   dfdx = @(x)[df1dx1(x),df1dx2(x);df2dx1(x),df2dx2(x)];
   x = x0;
   error = 1;
   while error >= e_d
    fx = f(x); dfx = dfdx(x);
    x = x - dfx\fx;
    error = norm(fx,Inf);
   end
   vin(i) = Vo*cos(t);
   vout(i) = V1+x(1);
end
plot(vout,vin)
grid on

[Ron]

Thanks for your help! This makes way more sense now