4.m

来自「利用龙格库塔法（等步长」· M 代码 · 共 47 行

clear all
a=0;b=10;
B=[0.1,0.5,1,10];
for j=1:4
    v=B(j);
    N=40
    h=(b-a)/(N+1)
    k=0.5;
    c=1;
    p=1;
    q0=50;
        T0=400;
    Tout=300;
    kv=10;
    l=1;m=3;
    xp=a:h:b;
    n=1;
    for x=a+h:h:l
        f(n)=0;
        n=n+1;
    end
    for x=l+h:h:m
     f(n)=q0*sin((x-l)/(m-l)*pi)*2*h^2;
    n=n+1;
    end
    for x=m+h:h:b
        f(n)=0
        n=n+1;
    end
  f(1)=(2*k+v*h*p*c)*T0;
    f(n-1)=(2*kv*h-v*p*c*kv*h^2*2/k)*Tout;
    A=zeros(N,N);
    A(1,1)=1;
    for i=2:N-1
        A(i,i-1)=-2*k-v*p*c*h;
        A(i,i)=4*k;
        A(i,i+1)=-2*k+v*p*c*h;
    end
    A(N,N-1)=-4*k;
    A(N,N)=4*k+4*h*kv-2*v*p*c*h^2/k
        u=A\f';
    up=[0,u',0]

    plot(xp,up)
    hold on
end