小弟求助各位大神，matlab解偏微分方程

小弟想求解这个偏微分方程组。

clear
clc
[FileName,PathName] = uigetfile('*.xls','Select the data code file');
FinalFileName=strcat(PathName,FileName);
Data=xlsread(FinalFileName);
global l yipuxilong EI Vc T0 rou D d ws ww H w Cd Tw Lw kesei k Cl mba lamdab1 lamdab2 lamdab3 lamdac1 lamdac2 lamdac3 C1 C2 C3 u Klz Kd x n1 n2 n3
n1=1;
n2=2;
n3=3;
l=Data(7,1);
x=1;
yipuxilong=Data(9,1);
EI=Data(1,1);
Vc=Data(11,1);
T0=Data(10,1);
rou=Data(13,1);
D=Data(14,1);
d=Data(22,1);
ws=Data(15,1);
ww=Data(20,1);
H=Data(19,1);
w=Data(18,1);
Cd=Data(17,1);
Tw=2*pi/ww;
Lw=Data(21,1);
kesei=Data(8,1);
k=2*pi/Lw;
x=1;
Cl=Data(16,1);
mba=Data(2,1)+Data(3,1);
lamdab1=(n1*pi/l).^4*(EI/mba);
lamdab2=(n2*pi/l).^4*(EI/mba);
lamdab3=(n3*pi/l).^4*(EI/mba);
lamdac1=(n1*pi/l).^4*(T0/mba);
lamdac2=(n2*pi/l).^4*(T0/mba)
lamdac3=(n3*pi/l).^4*(T0/mba)
C1=mba*sqrt((lamdab1+lamdac1))*yipuxilong;
C2=mba*sqrt((lamdab2+lamdac2))*yipuxilong;
C3=mba*sqrt((lamdab3+lamdac3))*yipuxilong;
% syms z t final y Y2 ideal
% zpie=z-(l+d);
% u=pi*H/Tw*exp(k*zpie)*cos(k*x+ww*t);
% Klz=0.5*rou*D*(Vc+u);
Kd=0.5*rou*D;
% calcu=2*Cl/(l*mba)*int(Klz*sin(n1*pi*z/l),z,0,l);
%
% Dn=(0.5*rou*D/2)*Cd*sign(Y2)*(Y2.^2)*(int(sin(n1*pi*z/l),z,0,l)).^3
%
% diff(Y2,t)=calcu*cos(ws*t)-(Cn/mba)*Y2-(lamdabn+lamdacn*(1+kesei*cos(w*t)))-2*Kd*Cd/(pi*mba)*Dn
%
%
% wn=EI*pi.^4./(mba*l.^4)+T0*pi./(mba*l.^2)
[T1,F1]=ode45('test_fun',[0 1],[0 0 0 0 0 0]);%测试时改变test_fun的函数维数，别忘记改变初始值的维数
plot(T1,F1)
function dy=test_fun(t,y)
dy = zeros(2,1);%初始化列向量
% dy(1) = y(1) - 2.*x./y(1)+y(1)+y(2)+y(3);
% dy(2) = y(2) - 2.*x./y(2)+y(1)+y(2)+y(3);
% dy(3) = y(3) - 2.*x./y(3)+y(1)+y(2)+y(3);
% global l yipuxilong EI Vc T0 rou D d ws ww H w Cd Tw Lw kesei k Cl mba lamdab1 lamdab2 lamdab3 lamdac1 lamdac2 lamdac3 C1 C2 C3 u Klz Kd x n1 n2 n3
% dy(1)=y(2);
% dy(2)=2*Cl./(l*mba)*cos(ws*t)*int(0.5*rou*D*(Vc+(pi*H./Tw*exp(k*(z-(l+d)))*cos(k*x+ww*t))).^2.*sin(n1.*pi*z./l),z,0,l)-(lamdab1+lamdac1*(1+kesei*cos(w*t)))*y(1)-C1./mba*y(2)+2*Kd*Cd./(pi*mba)*(Kd*Cd*int(sign(y(2).*sin(n1.*pi.*z./l)+y(4).*sin(n2.*pi.*z./l+y(6).*sin(n3.*pi.*z./l))).*(y(2).*sin(n1.*pi.*z./l)+y(4).*sin(n2.*pi.*z./l)+y(6).*sin(n3.*pi.*z./l)).^2.*sin(n1.*pi.*z./l),z,0,l));
% dy(3)=y(4);
% dy(4)=2*Cl./(l*mba)*cos(ws*t)*int(0.5*rou*D*(Vc+(pi*H./Tw*exp(k*(z-(l+d)))*cos(k*x+ww*t))).^2.*sin(n2.*pi*z./l),z,0,l)-(lamdab2+lamdac2*(1+kesei*cos(w*t)))*y(3)-C2./mba*y(4)+2*Kd*Cd./(pi*mba)*(Kd*Cd*int(sign(y(2).*sin(n1.*pi.*z./l)+y(4).*sin(n2.*pi.*z./l+y(6).*sin(n3.*pi.*z./l))).*(y(2).*sin(n1.*pi.*z./l)+y(4).*sin(n2.*pi.*z./l)+y(6).*sin(n3.*pi.*z./l)).^2.*sin(n2.*pi.*z./l),z,0,l));
% dy(5)=y(6);
% dy(6)=2*Cl./(l*mba)*cos(ws*t)*int(0.5*rou*D*(Vc+(pi*H./Tw*exp(k*(z-(l+d)))*cos(k*x+ww*t))).^2.*sin(n3.*pi*z./l),z,0,l)-(lamdab3+lamdac3*(1+kesei*cos(w*t)))*y(5)-C3./mba*y(6)+2*Kd*Cd./(pi*mba)*(Kd*Cd*int(sign(y(2).*sin(n1.*pi.*z./l)+y(4).*sin(n2.*pi.*z./l+y(6).*sin(n3.*pi.*z./l))).*(y(2).*sin(n1.*pi.*z./l)+y(4).*sin(n2.*pi.*z./l)+y(6).*sin(n3.*pi.*z./l)).^2.*sin(n3.*pi.*z./l),z,0,l));
%
global l yipuxilong EI Vc T0 rou D d ws ww H w Cd Tw Lw kesei k Cl mba lamdab1 lamdab2 lamdab3 lamdac1 lamdac2 lamdac3 C1 C2 C3 u Klz Kd x n1 n2 n3
syms z
dy(1)=y(2);
dy(2)=0.1304*Cl*cos(0.869*t)-0.00156*y(2)-(0.023+0.733*(1+0.328*cos(1.74*t)))*y(1)+0.1819*Cd*Kd*Cd*int(sign(y(2)*sin(pi*z./l)+y(4).*sin(n2*pi*z./l)+y(6).*sin(n3*pi*z./l)).*(y(2)*sin(pi*z./l)+y(4).*sin(n2*pi*z./l)+y(6).*sin(n3*pi*z./l)).^2*sin(n1*pi*z./l),z,0,l);
dy(3)=y(4);
dy(4)=-0.0416*Cl*cos(ws*t)-0.00326*y(4)-(0.3697+2.931*(1+kesei*cos(w*t)))*y(3)+0.1819*Cd*Kd*Cd*int(sign(y(2)*sin(pi*z./l)+y(4).*sin(n2*pi*z./l)+y(6).*sin(n3*pi*z./l)).*(y(2)*sin(pi*z./l)+y(4).*sin(n2*pi*z./l)+y(6).*sin(n3*pi*z./l)).^2*sin(n2*pi*z./l),z,0,l);
dy(5)=y(6);
dy(6)=0.0467*Cl*cos(ws*t)-0.00522*y(6)-(1.921+6.595*(1+kesei*cos(w*t)))*y(5)+0.1819*Cd*Kd*Cd*int(sign(y(2)*sin(pi*z./l)+y(4).*sin(n2*pi*z./l)+y(6).*sin(n3*pi*z./l)).*(y(2)*sin(pi*z./l)+y(4).*sin(n2*pi*z./l)+y(6).*sin(n3*pi*z./l)).^2*sin(n3*pi*z./l),z,0,l);
function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）
n=floor((b-a)/h);%求步数
x(1)=a;%时间起点
y(:,1)=y0;%赋初值，可以是向量，但是要注意维数
for ii=1:n
x(ii+1)=x(ii)+h;
k1=ufunc(x(ii),y(:,ii));
k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);
k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);
k4=ufunc(x(ii)+h,y(:,ii)+h*k3);
y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;
%按照龙格库塔方法进行数值求解
end
这是我参照网上编写的程序，求解不出来。请各位大神看一下问题出在哪里了呢？
结果应该是这样子的。