EH= 0.67;
EL= -0.33;

syms Va Vb Js n

T=1/100;
w=2*pi/T;
R23= -(R2+R3)/(R2*R3);
e1=Va+(K*(Vb-Va)*R3)/(R23)(R2+R3);
e2=(Vb-Va)/R23+Vb*i*n*w*C2-Js;
sol=solve(e1,e2,'Va,Vb');
Va=sol.Va
Vb=sol.Vb

iR23= (Vb-Va)/R23
iR2= -(iR23*R3)/(R2+R3)

HR23=iR23/Js;
HR2=-(H23*R3)/(R2+R3)


alfa = 0.05; %coefficiente di attenuazione

N=round(sqrt(500^2-(alfa^2)*1000^2)/(3*alfa*w));
clear mHR2
for n=1:N
    mHR2(n)=abs(subs(HR2));
end

figure, plot(mHR2,'o'), grid on, dock

clear n
syms n

%Serie di fourier, costruire l'ingresso.
t1=0.33*T;

I0 = 0; %Valore medio nullo
An_sym=2/T*(int(EH*cos(n*w*t),t,0,t1)+int(EL*cos(n*w*t),t,t1,T));
Bn_sym=2/T*(int(EH*sin(n*w*t),t,0,t1)+int(EL*sin(n*w*t),t,t1,T));

J_t=I0; %termine costante ingresso
iR2_0= subs(HR2,'n',0)*E0; %termine costante uscita
iR2_t= iR2_0;

    vet_FCn= zeros(N,1);
    vet_HR2=zeros(N,1);
    vet_FiR2=zeros(N,1);
    vet_N=[1:1:N]';
for n=1:N
    %ingresso
    An=subs(An_sym); %l'integrale è in n, con la subs ad ogni iterazione sostituisce la n.
    Bn=subs(Bn_sym);
    FCn=An-i*Bn;
    J_t=J_t+vpa(abs(FCn)*cos(w*n*t+angle(FCn))); %serie fourier ingresso

    %uscita iR2
    HR2n=subs(HR2);
    FiR2_Cn=HR2n* FCn;
    iR2_n= vpa(abs(FiR2_Cn)*cos(n*w*t+angle(FiR2_Cn))); %componenti uscita

    iR2_t= iR2_t + iR2_n; %serie fourier uscita

    vet_FCn(n)= abs(FCn);
    vet_HR2(n)=abs(HR2n);
    vet_FiR2(n)=abs(FiR2_Cn);
end

figure (4), hj= ezplot(J_t,[0,2*T]),grid on, dock, axis auto, hold on
figure (4), hr= ezplot(iR2_t,[0,2*T]),grid on, dock, axis auto

set(hj,'color','r'), set(hr,'color','b')