MATLAB: Can anyone help me rectifying the error shown while executing this program

                kermat(:,col) = Kvals(:,col)*wt(col);

                [told,solnold] = mod_simp(n);

                [soln,errest,cond] = 
                Fie(lambda,a,b,behavior,@kernel,@RHS,AbsTol,RelTol);

                TestSmooth(1,0,1)
       %==================The program is as follows;
       function TestAll
       disp('SMOOTH KERNEL')
       TestSmooth(1,0,1)
       disp(' ')
       disp('********************************************************************')
       disp(' ')
       end
       %==========================================================================

       function TestSmooth(lambda,a,b,AbsTol,RelTol)
       behavior =1;
       if nargin < 5 isempty(AbsTol)
       AbsTol = 1e-6;
       end
       if nargin < 6 isempty(RelTol)
       RelTol = 1e-3;
       end
       % Compute and plot the solution:
       [soln,errest,cond] = Fie(lambda,a,b,behavior,@kernel,@RHS,AbsTol,RelTol);
       t = soln.s; sol = soln.x;
       nfinal = length(t) - 1;
       % Interpolate the solution for assessing the error and if necessary,
       % use to get a smooth graph.
       tint = linspace(a,b,150);
       xint = ntrpFie(soln,tint);
       if nfinal < 150
       plot(tint,xint)
       else
       plot(t,sol)
       end
       % Study the error and condition of the integral equation.
       true_error = norm(true_soln(t) - sol,inf);
       max_error = norm(true_soln(tint)-xint,inf);
       disp(' ')
       disp(['Condition number = ',num2str(cond)])
       disp(' ')
       fprintf('Approximate bound on error at nodes = %6.1e\n',errest)
       fprintf('Actual error at nodes = %6.1e\n',true_error)
       fprintf('\nActual error at 150 equally spaced points = %6.1e\n',max_error)
       function kst = kernel(s,t)
       nn=16;
       height=20e-6; % half-height of the plate
       aa=10e-6;
       F3=@(ss)(1-exp(-2.*(ss./aa).*height))./((1+exp(-2.*
       (ss./aa).*height)));%F3(s/aa)
       sdata=0:0.1:1;tdata=0:0.1:1;
       for ii=1:length(sdata);
       s=sdata(ii);
       for jj=1:length(tdata);
       t=tdata(jj);
       f21=@(ss)ss.*(F3(ss)-1).*besselj(0,ss.*s.^2).*besselj(0,ss.*t.^2);
       f2=@(ss)-2.*t.*s.*t.*f21(ss);
       I(ii,jj)=galag(f2,nn);
       end
       end
       J=sum(I); % kernel
       kst=sum(J);
       function ts = true_soln(t)
       ts = 0;
       end
       function rs = RHS(s)
       rs= s ;
       end
       end
       function s = galag(func,n)
       if (n==2)|(n==4)|(n==8)|n(n==16)
       c=zeros(16,3); t=zeros(16,3);
       c(1:2,1)=[1.533326033   4.45095733]; % Cs are the weights
       c(1:4,2)=[0.832739124   2.048102438 3.631146306 6.487145084];
       c(1:8,3)=[0.43772341    1.033869348 1.669709766 2.376924702 3.208540913 
       4.268575511 5.818083369 8.906226215];
       c(:,4)=[0.225036315 0.525836053 0.831961392 1.146099241 1.471751317 
       1.813134687 2.17551752 2.56576275   2.993215086 3.471234483 4.020044086 
       4.672516608 5.487420658 6.585361233 8.276357984 11.82427755];
       t(1:2,1)=[0.585786438   3.414213562]; % ts are the nodes
       t(1:4,2)=[0.32254769    1.745761101 4.536620297 9.395070912];
       t(1:8,3)=[0.170279632   0.903701777 2.25108663  4.26670017  7.045905402 
       10.75851601 15.74067864 22.86313174];
       t(:,4)=[0.08764941  0.462696329 1.141057775 2.129283645 3.437086634 
       5.078018615 7.070338535 9.438314336 12.21422337 15.44152737 19.18015686 
       23.51590569 28.57872974 34.5833987  41.94045265 51.70116034];
       j=1;
        while j<=3
        if 2^j==n; break
        else
        j=j+1;
        end
        end
        s=0;
        for k=1:n
        x=t(k,j);y=feval(func,x);
        s=s+c(k,j)*y;
        end
        end
        end
        end
        %==========================================================================
        function [sol,errest,cond] = 
        Fie(lambda,a,b,behavior,kernel,RHS,AbsTol,RelTol)
        % Create and initialize solution structure sol.
        n = 2^3;
        [told,solnold] = mod_simp(n);
        if behavior == 0
        sol.s = (1 - told) ./ told;
        else
        sol.s = told;
        end sol.x = solnold;
        sol.lambda = lambda;
        sol.kernel = kernel;
        sol.behavior = behavior;
        if behavior == 4
        sol.alpha = alpha;
        end
        sol.RHS = RHS;
        max_m = 9;
        for m = 4:max_m
        n = 2^m;
        [t,soln] = mod_simp(n);   % Solve with n subdivisions of [a,b].
        % Calculate the norm of the difference between the current numerical
        % solution and the preceding one. In several cases we must interpolate
        % to get approximations on the same mesh.
        delta = norm(soln(1:2:n+1) - solnold,inf);
        % Update solution structure.
        sol.s = t;
        sol.x = soln;
        if m > 4
        % Estimate contraction rate and use to predict error.  
        % Be conservative and do not use if not contracting.
        rate = min(0.5,max(delta/deltaold,0.0625));
        end
        errest = (rate/(1 - rate))*delta; 
        if errest < (atol + rtol*norm(soln,inf))
        break;
        end  
        deltaold = delta;
        if m == max_m
        warning('Fie:Failure','Failed to pass the error test.')
        end
        end
        % Inexpensive approximation of condition of matrix in 1-norm.
        if nargout == 3
        cond = condest(kermat);
        end
        function [t,soln] = mod_simp(n)
        % Solve with n subdivisions of [a,b] using a quadrature rule. Simpson's
        % rule is used for smooth kernels.
        h = (b - a)/n;
        t = linspace(a,b,n+1)';
        % Weights for Simpson's rule.  Assumes n is even.
        wt = ones(1,n+1); wt(3:2:n-1) = 2; wt(2:2:n) = 4;
        wt = (h/3)*wt;
        [S,T] = meshgrid(t,t);
        rhs = RHS;
        Kvals = kernel(S,T)';
        Msize = n + 1;
        kermat = zeros(Msize);
        for col = 1:n+1
            kermat(:,col) = Kvals(:,col)*wt(col);
        end
        for j = 1:Msize
        kermat(j,j) = kermat(j,j) - lambda;
        end
        soln = kermat \ (-rhs);
        end % mod_simp
        end