lab61.m

利用龙格库塔法（等步长

figure(1)
clear all

global N hx

N=40;
hx=1/N;
x=0:hx:1;
M=200;
ht=5/M;
t=0:ht:5;
c1=ht^2/(2*hx^2)+ht/(2*hx)-(0.3*ht^2)/(2*hx);
c2=1-ht^2/hx^2-0.3*ht+0.3^2*ht^2/2;
c3=ht^2/(2*hx^2)-ht/(2*hx)+(0.3*ht^2)/(2*hx);
c4=0.3*ht*5/3-0.3^2*5/3*ht^2/2;
A=zeros(M+2,M+2);
for j=1:N+1
    A(j,1)=5;
end
for i=1:M+1
    if t(i)<=1/6&&t(i)>=0 
        A(1,i)=(5+95*sin(pi*3*t(i)));
    elseif t(i)<=4/3&&t(i)>=1/6
        A(1,i)=100;
    else
        A(1,i)=(100+5*sin(pi*(3*t(i)-4)));
    end
end

for k=1:M
    for j=2:N+1
       A(j,k+1)=c1*A(j-1,k)+c2*A(j,k)+c3*A(j+1,k)+c4;
       A(j,M+2)=A(j,M+1);
    end
end
for k=1:M+1
    for j=1:N+1
        B(k,j)=A(j,k);
    end
end


mesh(x,t,B)
xlabel('rod length variable x')
ylabel('time variable t')
zlabel('temperature variable T')