2.m

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

clear all
a=1;b=3;
B=[0.1,0.5,1,10];
for j=1:4
    v=B(j);
    N=40;
    h=(b-a)/N;
    k=0.5;
    c=1;
    p=1;
    q0=50;
    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)
    n=n+1;
    end
    for x=m+h:h:b-h
        f(n)=0
        n=n+1;
    end
    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-2)=-k;
    A(N,N-1)=-2*h*k*v;
    A(N,N)=-k
        u=A\f';
    up=[0,u',0]

    plot(xp,up)
    hold on
end