I am trying to solve a set of ODEs using the Harmonic Balance method. In order to do this, I need to compute the Jacobian of the set of equations. However I am very confused regarding the dimensions of this Jacobian matrix. In principle, given that I have 3 equations, the Jacobian would be a $3\times 3$ matrix. However, in Harmonic Balance, the optimisation is done over the entire set of frequencies at any given iteration. This means that in practice, the matrix containing the ODEs, instead of being $3\times 1$, it has $3\times N_\omega$, where $N_\omega$ is the number of frequency coeficients that I am considering. Does this means that I have to compute $N_\omega$ Jacobians at every iteration, one per frequency? I leave below the MATLAB code I have devoleped so far, without the calculation of the Jacobian, in order to ilustrate my problem. How can I compute this Jacobian?
fs = 32; w = 2*pi*fs; N=10*64; X=zeros(N,1); X(fs)=N; X(1)=20*N; X = X'; s = rand(3,length(X)); tol = 1e-3; max_iter = 50; i=0; e = 9; while(e>tol && i<max_iter) e = 1i*w*s-F(s,X); J = jacobian(F,s) % <-Here is the doubt, how do I compute this Jacobian s = s - (1i*w - J(s, Gm, Gds))\e; i = i + 1; end function y = F(s, X) VDD = zeros(1,length(X)); VDD(1) = 10; Ld = 120e-6; Cg = 50e-12; Cb = 40e-12; RS = 50; RL = 6.6; vG = s(1,:); vD = s(2,:); iD = fft(10*(1/2.*(ifft(vG)-3)+1/20.*log(2*cosh(10*(ifft(vG)-3)))).*(1+0.003.*ifft(vD)).*tanh(ifft(vD))); y = [ 1/Cg*((X-s(1,:))/RS+iD+s(3,:)+s(2,:)/RL); 1/Cg*(X-s(1,:))/RS-(Cb+Cg)/(Cb*Cg)*(iD+s(3,:)+s(2,:)/RL); (s(2,:)-VDD)/Ld; ]; end