lab6.m

利用龙格库塔法（等步长

figure(2)
clear all

global N hx
N=40;
hx=1/N;
x=hx:hx:1;
u0=zeros(1,N)+5;


t0=0;tend=5;
time=[t0];
u=u0';
M=420;
t=t0;
result=u0;
for k=1:M
    ht=(tend-t0)/M;
    u=u+ht*fpar(t,u);
    t=t+ht;
    result=[result;u'];
    time=[time;t];
end
for i=1:size(time)
    if time(i)<=1/6&&time(i)>=0 
        A(i)=(5+95*sin(pi*3*time(i)));
    elseif time(i)<=4/3&&time(i)>=1/6
        A(i)=100;
    else
        A(i)=(100+5*sin(pi*(3*time(i)-4)));
    end
end
resultp=[A' result];
xp=0:hx:1;
mesh(xp,time,resultp)
xlabel('rod length variable x')
ylabel('time variable t')
zlabel('temperature variable T')