Исправить ошибку решение диф. уравнения

%  file  ka1_e.m
 
clc
clf
a0=100;
s1='                     Введите      ';
s2='                   значение l    ';
disp(s1),disp(s2)
l=input('l = ');
s3='     параметры  входного гармонического воздействия   ';
disp(s1),disp(s3)
u0=input('u0 =');
w0=input('w0 =');
s4='    параметры  шумового гармонического воздействия     ';
disp(s1),disp(s4)
u1=input(' u1 = ');
w1=input(' w1 = ');
s5='           интервал  времени : tm-tm1';
disp(s1),disp(s5)
tm=input('tm =');
tm1=input('tm1 =');
 
t0=tm;
tfinal=tm1;
y0=[0 0 0 0 0 0]'; %  начальные значения вектора координат системы, см.
% файл file  vdp1.m
tol=0.00001; %  шаг изменения временной координаты – точность % представления данных
[t,y]=ode23('vdp1',t0,tfinal,y0,tol);
clf
 
subplot(121)
plot(t,y(:,1),'-g')
grid
title('   Координата  Xd   ')
xlabel(' t  ,  с')
ylabel(' xd')
subplot(122)
plot(t,y(:,2),'-w')
grid
title('    Координата  Xm  ')
xlabel(' t ,  с')
ylabel(' xm')
pause
clf
 
k00=(-2*l)*y(:,6);        %  бв®«ЎҐж
k0=k00';                  %  бва®Є 
u=u0*sin(w0*t);           %  бв®«ЎҐж
u00=u';                   %  бва®Є 
u2=u1*sin(w1*t);          %  бв®«ЎҐж
u22=u2';                  %  бва®Є 
yd=y(:,1)+u2;          %  бв®«ЎҐж
ydd=yd';                  %  бва®Є 
 
n=length(k0);             %  ¤«Ё¬  ўҐЄв®а  k0
for i=1:n
a(i)=100;                 %  бва®Є 
end
a00=a';                   %  бв®«ЎҐж
y2=y(:,2)';               %  бва®Є 
um=(ydd+y2.*k0).*((a-k0)./a);    %  бва®Є 
umm=um';                         %  бв®«ЎҐж
 
subplot(121)
plot(t,yd,'-g')
grid
title(' Координата  ( Yd ) ')
xlabel('  t,  б')
ylabel(' yd')
subplot(122)
plot(t,umm,'-r')
grid
title(' Координата ( Um )')
xlabel('  t,  б')
ylabel(' um')
pause
clf
 
subplot(221)
plot(t,y(:,3),'-g')
title('  Координата   u10')
xlabel(' t , б')
ylabel(' u10')
grid
subplot(222)
plot(t,y(:,4),'-g')
title('  Координата   z1')
xlabel(' t ,б')
ylabel(' z1')
grid
subplot(223)
plot(t,y(:,5),'-g')
title('  Координата   q10')
xlabel(' t , с')
ylabel(' q10')
grid
subplot(224)
plot(t,y(:,6),'-g')
grid
title('  Координата  w')
xlabel(' t ,  с')
ylabel(' w')
grid
pause
clf
 
plot(t,k00,'-g')
grid
title(' коэффициент  подстройки ')
xlabel(' t , с')
ylabel(' k0')
pause
clf
 
ei=u00-um;                    %   бва®Є 
eii=ei';                      %   бв®«ЎҐж
y4=y(:,4)';                   %   бва®Є 
e1=(ydd-(a-k0).*y2+k0.*y4).*((a-k0)./a);   % бва®Є 
e11=e1';                      %  бв®«ЎҐж
 
subplot(121)
plot(t,eii,'-g')
grid
title('  Погрешность Е ист ')
xlabel(' t  , с')
ylabel('  e ист')
subplot(122)
plot(t,e11,'-g')
grid
title('  Погрешность Е ип')
xlabel('  t  , с')
ylabel(' e ип')
pause
clf

% file  vdp1.m
function yprime=vdp1(t,y);
a0=100;
 
l=100000000;
 
u0=1;w0=100;
u1=0.05;
w1=1000;
 
u=u0*sin(w0*t);
u2=u1*sin(w1*t);
yprime(1)=u-a0*y(1);                   % xd
yd=y(1)+u2;
k0=(-2*l)*y(6);
yprime(2)=yd+(k0-a0)*y(2);             % xm
ym=(a0-k0)*y(2);
yprime(3)=y(2)+(k0-a0)*y(3);           % u10
e=yd-ym;
yprime(4)=e+(k0-a0)*y(4);              % z1
yprime(5)=y(4)-y(3)+(k0-a0)*y(5);   % q10
e1=(e+k0*y(4))*((a0-k0)/a0);
yprime(6)=e1*(y(4)+k0*y(5)-y(3));   % w
yprime=yprime';