rk4.txt

龙哥库塔4阶函数的计算

function R=rk4(f,t,y0,h)

%Input    - f is the function 
%            - t is the right endpoint
%            - yO is the initial condition y(O)
%            - h is the index number
%Output - R = Y(M+1) 
%              and Y(M+1) is the solution 


M=t/h;
T=zeros(1,M+1);
Y=zeros(1,M+1);
T=0:h:t;
Y(1)=y0;
for j=1:M
   k1=h*f(T(j),Y(j));
   k2=h*f(T(j)+h/2,Y(j)+k1/2);
   k3=h*f(T(j)+h/2,Y(j)+k2/2);
   k4=h*f(T(j)+h,Y(j)+k3);
   Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;
end

R=Y(M+1);

return