用MATLAB解偏微分方程

set(0,'defaultfigurecolor','w')
h=0.1;
N=30;
dt=0.0001;
M=10000;
A=dt/(h^2);
U=zeros([N+1,M+1]);
Space=0:h:(N)*h;


for k=1:M+1
    U(1,k) = 0.0;
    U(N+1,k) = 0.0;
end


for i =1:N
    U(i,1)=4*(i-1)*h*(3-(i-1)*h);
end


for k=1:M
    for i=2:N
        U(i,k+1)=A*U(i+1,k)+(1-2*A)*U(i,k)+A*U(i-1,k);
    end
end
figure;
plot(Space,U(:,1), 'g-')
hold on;
plot(Space,U(:,3000), 'b-')
hold on;
plot(Space,U(:,6000), 'k-')
hold on;
plot(Space,U(:,9000),'r-')
hold on;
plot(Space,U(:,10000),'y-')
ylabel('u(x,t)')
xlabel('x')
xlim([0,3])
ylim([-2,10])
legend('t=0','t=0.3','t=0.6','t=0.9','t=1')
figure;
contourf(U)
colorbar