Ошибка при присвоении элементу матрицы вычисленного значения

X=hw01_int_digit(t,y,'левый');

function [ output_args ] = Untitled2( input_args )
%основная функция
clc;
a=0.1; b=0.7; eps=1e-6;
n=1000
 
 h=(b-a)/n;
 t=a:h:b;
 y=zeros(1,n+1); 
  for j=1:n+1
     y(j)=hw01_func(t(j));
  end  
  S1=zeros(1,4);
  S1(1,1)=n; S1(1,2)=hw01_int_digit(t,y,'левый');
  S1(1,3)=hw01_int_digit(t,y,'правый');
  S1(1,4)= hw01_int_digit(t,y,'трапеций');
 
 S1
 
 function y= hw01_func(t)
  y= 0.89*t + 0.13* sin(t)-1.16*t^2; 
  
 function f=hw01_int_digit(t,y,s) 
% t -массив времен, y-массив знач, s назв метода интегрирования
 f=0; n=length(t);
switch s
case 'левый'
    for i=1:n-1
        f=f+(t(i+1)-t(i))*y(i);
     end;
case 'правый'
    for i=1:n-1
        f=f+(t(i+1)-t(i))*y(i+1);
     end;
case 'трапеций'
     for i=1:n-1
       f=f+(t(i+1)-t(i))*(y(i)+y(i+1))*0.5;   
     end;
end

%основная функция
clc;
a=0.1;
b=0.7;
eps=1e-6;
n=1000 
t=linspace(a,b,n);
y=0.89*t+0.13*sin(t)-1.16*t.^2; 
S1=zeros(1,4);
S1(1,1)=n;
S1(1,2)=hw01_int_digit(t,y,'левый');
S1(1,3)=hw01_int_digit(t,y,'правый');
S1(1,4)= hw01_int_digit(t,y,'трапеций');
S1

function f=hw01_int_digit(T,Y,s) 
% T - массив времени, Y-массив знач, s - назв метода интегрирования
n=length(T);
TT=circshift(T',-1)-T';
switch s
  case 'левый'
    YY=Y(1:end-1);
  case 'правый'
    YY=Y(2:end);
  case 'трапеций'
    YY=(Y'+circshift(Y',-1))'/2;
    YY=YY(1:end-1);
end
f=sum(TT(1:end-1).*YY');
end