# Estimating the Jacobian in Harmonic Balance Method

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