The following code is to solve a hydrodynamic stability problem. It is a shooting method implemented in Octave software

Question

I have the following code written in Octave. It is supposed to solve a coupled differential equation (DE) dealing with hydrodynamic stability analysis. The code works fine without imaginary numbers and gives excellent result. However, as most DE in stability analysis contain imaginary numbers, I want my code to work when imaginary numbers are included in the DE.

function r = f1 (y,k)
NN = 2000;
Pr = 0.7;
eta = 0;
Ra = y(1); xi = y(2); zeta = y(3);
f = @(x, t) [x(2);
x(3);
x(4);
(2*k^2+eta/Pr)*x(3)-(k^2+eta/Pr)*k^2*x(1)+Ra*k^2*x(5);
x(6);
(eta+k^2)*x(5)-i*x(1)];
t = linspace (0, 1, NN);
lsode_options("absolute tolerance", 1e-16);
lsode_options("relative tolerance", 1e-15);
x = lsode (f, [0; 0; xi; zeta; 0; 1], t);
r(1) = x(NN,1);
r(2) = x(NN,2);
r(3) = x(NN,5);
endfunction


fileID = fopen('file1.txt','w');
kk = linspace(kmin,kmax,Nmax);
yy = [1708; 1; -1];
yyold = yy;
options=optimset('MaxIter',1e5,'TolFun',1e-16,'ComplexEqn','on');
hold on;
for ll = 1:Nmax
k = kk(ll);
yy = fsolve(@(y) f1(y,k),yyold,options);
yyold = yy;
format long;
printf("k = %.15f\n", k);
printf("Ra = %.15f\n", yy(1));
plot(k,yy(1))
fflush(stdout);
fprintf(fileID,'( %.10f , %.10f )\n', k, yy(1));
endfor
fclose(fileID);

In the optimset I have put the option of ComplexEqn. But to no avail. I want to solve the given equations seamlessly.