matlab:使用四阶龙格库塔方法求解微分方程

%书籍：常用数值算法及其matlab实现
%第10章 常微分方程初值问题的数值解法
%四阶龙格库塔方法
function S = rungek4(fun, x0, xn, y0, h)
%fun:微分方程的右表达式
%x0, xn为区间
%y0 为初值
M = floor(xn-x0)/h ;      %离散点的个数M+1
T =zeros(1, M+1); Y =zeros(1, M+1);         %行向量
T(1) = x0;
Y(1) = y0;
for i = 1:M
    K1 =  feval(fun, T(i) ,Y(i));
    K2 =  feval(fun, T(i)+1/2*h ,Y(i)+1/2* h*K1);
    K3 =  feval(fun, T(i)+1/2*h ,Y(i)+1/2* h*K2);
    K4 =  feval(fun, T(i)+ h ,Y(i)+ h*K3);   
    Y(i+1) = Y(i) +h/6 *(K1 + 2*K2 + 2*K3 + K4);
    T(i+1) = T(i) + h;
end
S = [T' Y'];            %E是M+1行，2列
主函数
%书籍：常用数值算法及其matlab实现
%第10章 常微分方程初值问题的数值解法
%四阶龙格库塔，例10.5
%function S = rungek4(fun, x0, xn, y0, h)
clear all;clc;close  all;
%fun =@(x,y)2*(9*x^2 - x + 9)/(1+ x^2)^2 - 18*y;
fun = inline('2*(9*x^2 - x + 9)/(1+ x^2)^2 - 18*y');
a = 0; b =1;