try running here, will give errors: https://rextester.com/l/octave_online_compiler

LHTSS.m

% data taken from Excel sheet described in:
% http://scottburns.us/reflectance-curves-from-srgb-10/
T = [
0.0000646936115727633 0.000219415369171578 0.00112060228414359 0.00376670730427686 0.0118808497572766 0.0232870228938867 0.0345602796797156 0.0372247180152918 0.0324191842208867 0.0212337349018611 0.0104912522835777 0.00329591973705558 0.000507047802540891 0.000948697853868474 0.00627387448845597 0.0168650445840847 0.0286903641895679 0.0426758762490725 0.0562561504260008 0.0694721289967602 0.0830552220141023 0.0861282432155783 0.0904683927868683 0.0850059839999687 0.0709084366392777 0.0506301536932269 0.0354748461653679 0.0214687454102844 0.0125167687669176 0.00680475126078526 0.00346465215790157 0.00149764708248624 0.000769719667700118 0.000407378212832335 0.000169014616182123 0.0000952268887534793;
0.00000184433541764457 0.0000062054782702308 0.0000310103776744139 0.000104750996050908 0.000353649345357243 0.000951495123526191 0.00228232006613489 0.00420743392201395 0.00668896510747318 0.00988864251316196 0.015249831581587 0.0214188448516808 0.0334237633103485 0.0513112925264347 0.0704038388936896 0.0878408968669549 0.0942514030194481 0.0979591120948518 0.094154532672617 0.0867831869897857 0.078858499565938 0.0635282861874625 0.0537427564004085 0.0426471274206905 0.0316181374233466 0.0208857265390802 0.0138604556350511 0.00810284218307029 0.00463021767605804 0.002491442109212 0.00125933475912608 0.000541660024106255 0.000277959820700288 0.000147111734433903 0.0000610342686915558 0.0000343881801451621;
0.000305024750978023 0.00103683251144092 0.00531326877604233 0.0179548401495523 0.057079004340659 0.11365445199637 0.173363047597462 0.196211466514214 0.186087009289904 0.139953964010199 0.0891767523322851 0.0478974052884572 0.0281463269981882 0.0161380645679562 0.00775929533717298 0.00429625546625385 0.00200555920471153 0.000861492584272158 0.000369047917008248 0.000191433500712763 0.000149559313956664 0.0000923132295986905 0.0000681366166724671 0.0000288270841412222 0.0000157675750930075 0.00000394070233244055 0.00000158405207257727 0 0 0 0 0 0 0 0 0;
]
sRGB = [250.0, 92.0, 8.0]

% This is the Least Hyperbolic Tangent Slope Squared (LHTSS) algorithm for
% generating a "reasonable" reflectance curve from a given sRGB color triplet.
% The reflectance spans the wavelength range 380-730 nm in 10 nm increments.

% Solves min sum(z_i+1 - z_i)^2 s.t. T ((tanh(z)+1)/2) = rgb,
% using Lagrangian formulation and Newton's method.
% Reflectance will always be in the open interval (0->1).

% T is 3x36 matrix converting reflectance to D65-weighted linear rgb,
% sRGB is a 3 element vector of target D65 referenced sRGB values in 0-255 range,
% rho is a 36x1 vector of reconstructed reflectance values, all (0->1),

% Written by Scott Allen Burns, May 2019.
% Licensed under a Creative Commons Attribution-ShareAlike 4.0 International
% License (http://creativecommons.org/licenses/by-sa/4.0/).
% For more information, see http://scottburns.us/reflectance-curves-from-srgb/

% initialize outputs to zeros
rho=zeros(36,1);

% handle special case of (0,0,0)
if all(sRGB==0)
    rho=0.0001*ones(36,1);
    return
end

% handle special case of (255,255,255)
if all(sRGB==255)
    rho=ones(36,1);
    return
end

% 36x36 difference matrix for Jacobian
% having 4 on main diagonal and -2 on off diagonals,
% except first and last main diagonal are 2.
D=full(gallery('tridiag',36,-2,4,-2));
D(1,1)=2;
D(36,36)=2;

% compute target linear rgb values
sRGB=sRGB(:)/255; % convert to 0-1
rgb=zeros(3,1);
% remove gamma correction to get linear rgb
for i=1:3
    if sRGB(i)<0.04045
        rgb(i)=sRGB(i)/12.92;
    else
        rgb(i)=((sRGB(i)+0.055)/1.055)^2.4;
    end
end

% initialize
z=zeros(36,1); % starting point all rho=1/2
lambda=zeros(3,1); % starting Lagrange mult
maxit=100; % max number of iterations
ftol=1.0e-8; % function solution tolerance
count=0; % iteration counter

% Newton's method iteration
while count <= maxit
    d0 = (tanh(z) + 1)/2;
    d1 = diag((sech(z).^2)/2);
    d2 = diag(-sech(z).^2.*tanh(z));
    F = [D*z + d1*T'*lambda; T*d0 - rgb]; % 39x1 F vector
    J = [D + diag(d2*T'*lambda), d1*T'; T*d1, zeros(3)]; % 39x39 J matrix
    delta=J\(-F); % solve Newton system of equations J*delta = -F
    z=z+delta(1:36); % update z
    lambda=lambda+delta(37:39); % update lambda
    if all(abs(F)<ftol)
        % solution found
        rho=(tanh(z)+1)/2;
        return
    end
    count=count+1;
end
disp(['No solution found in ',num2str(maxit),' iterations.'])