运行这段在MATLAB中写的迎风法解双曲型PDE时稳定闪退

h=0.01;
r=0.5;
tao=r*h;
x=0:h:5;
t=0:tao:10;
length_t=length(t);
length_x=length(x);
U=zeros(length_t,length_x);
length_in=length_x-1;
A=diag((1-r)*ones(length_in,1))+r*diag(ones(length_in-1,1),-1);
[X,T]=meshgrid(x,t);
F=2*T.*sin(X)+T.^2.*cos(X);
for i=2:length(t)
    U(i,2:end)=A*U(i-1,2:end)'+tao*F(i,2:end)';
end
U_acc=T.^2.*sin(X);
figure(1)
subplot(1,2,1)
mesh(x,t,U)
title("解曲面")
subplot(1,2,2)
mesh(x,t,U_acc-U)
title("误差曲面")

发现是FOR循环内的矩阵转置处出现了问题，但是无法清晰定位具体发生了什么导致的问题，希望开发组能解决。