Error "not enough input arguments"

% Simulate supercontinuum generation for parameters similar
% to Fig.3 of Dudley et. al, RMP 78 1135 (2006)
% Written by J.C. Travers, M.H Frosz and J.M. Dudley (2009)
% Please cite this chapter in any publication using this code.
% Updates to this code are available at www.scgbook.info
n = 2^13;                   % number of grid points
twidth = 12.5;              % width of time window [ps]
c = 299792458*1e9/1e12;     % speed of light [nm/ps]
wavelength = 835;           % reference wavelength [nm]
w0 = (2.0*pi*c)/wavelength; % reference frequency [2*pi*THz]
T = linspace(-twidth/2, twidth/2, n); % time grid
% === input pulse
power = 10000;              % peak power of input [W]
t0 = 0.0284;                % duration of input [ps]
A = sqrt(power)*sech(T/t0); % input field [W^(1/2)]
% === fibre parameters
flength = 0.15;             % fibre length [m]
% betas = [beta2, beta3, ...] in units [ps^2/m, ps^3/m ...]
betas = [-11.830e-3, 8.1038e-5, -9.5205e-8, 2.0737e-10, ...
         -5.3943e-13, 1.3486e-15, -2.5495e-18, 3.0524e-21, ...
         -1.7140e-24];
gamma = 0.11;               % nonlinear coefficient [1/W/m]
loss = 0;                   % loss [dB/m]
% === Raman response
fr = 0.18;                  % fractional Raman contribution
tau1 = 0.0122; tau2 = 0.032;
RT = (tau1^2+tau2^2)/tau1/tau2^2*exp(-T/tau2).*sin(T/tau1);
RT(T<0) = 0;          % heaviside step function
%RT = RT/trapz(T,RT);  % normalise RT to unit integral
% === simulation parameters
nsaves = 200;     % number of length steps to save field at
% propagate field
[Z, AT, AW, W] = gnlse(T, A, w0, gamma, betas, loss, ...
                       fr, RT, flength, nsaves);
% === plot output
figure();
lIW = 10*log10(abs(AW).^2); % log scale spectral intensity
mlIW = max(max(lIW));       % max value, for scaling plot
WL = 2*pi*c./W; iis = (WL>400 & WL<1350); % wavelength grid
subplot(1,2,1);             
pcolor(WL(iis), Z, lIW(:,iis)); % plot as pseudocolor map
caxis([mlIW-40.0, mlIW]);  xlim([400,1350]); shading interp; 
xlabel('Wavelength / nm'); ylabel('Distance / m');

lIT = 10*log10(abs(AT).^2); % log scale temporal intensity
mlIT = max(max(lIT));       % max value, for scaling plot
subplot(1,2,2);
pcolor(T, Z, lIT);          % plot as pseudocolor map
caxis([mlIT-40.0, mlIT]);  xlim([-0.5,5]); shading interp;
xlabel('Delay / ps'); ylabel('Distance / m');