所解方程为  y’’= -y

主程序 M文件 

a=0;

b=8*pi;

n=400;

h=(b-a)/n;

xa=3;

ya=0;

T=a:h:b;

X=zeros(1,n+1);

Y=zeros(1,n+1);

X(1)=xa;

Y(1)=ya;

for i=1:n

    f1=f(T(i),X(i),Y(i));

    g1=g(T(i),X(i),Y(i));

    f2=f(T(i)+h/2,X(i)+f1*h/2,Y(i)+g1*h/2);

    g2=g(T(i)+h/2,X(i)+f1*h/2,Y(i)+g1*h/2);

    f3=f(T(i)+h/2,X(i)+f2*h/2,Y(i)+g2*h/2);

    g3=g(T(i)+h/2,X(i)+f2*h/2,Y(i)+g2*h/2);

    f4=f(T(i)+h,X(i)+f3*h,Y(i)+g3*h);

    g4=g(T(i)+h,X(i)+f3*h,Y(i)+g3*h);

    X(i+1)=X(i)+(f1+2*f2+2*f3+f4)*h/6;

    Y(i+1)=Y(i)+(g1+2*g2+2*g3+g4)*h/6;%四阶龙格-库塔方法

end

R1=[T',X',Y']

函数M文件

f.m

function y=f(x,y,z)

y=z; 

g.m

function y=g(x,y,z)

y=-y;

将两段M文件的程序分别放入M文件，即可运行。