matlab：二阶龙格库塔求解欧拉方程

%书籍：常用数值算法及其matlab实现%第10章 常微分方程初值问题的数值解法%二阶龙格库塔方法function S = heunsec(fun, x0, xn, y0, h)%fun:微分方程的右表达式%x0, xn为区间%y0 为初值M = floor(xn-x0)/h ;%离散点的个数M+1T =zeros(1, M+1); Y =zeros(1, M+1); %行向量T(1) = x0;Y(1) = y0;for i = 1:MK1 = feval(fun, T(i) ,Y(i));K2 = feval(fun, T(i)+2/3*h ,Y(i)+2/3* h*K1);Y(i+1) = Y(i) +h/4 *(K1 + 3*K2);T(i+1) = T(i) +h;endS = [T' Y']; %E是M+1行，2列主函数：%书籍：常用数值算法及其matlab实现%第10章 常微分方程初值问题的数值解法%二阶龙格库塔，例10.3%function S = heunsec(fun, x0, xn, y0, h)clear all;clc;close all;%fun =@(x,y)sin(x*y);fun = inline('sin(x*y)');a = -1; b =1;y0 = 0.5;h = 0.2S = heunsec(fun, a, b, y0, h)plot(S(:,1),S(:,2) ,'r*-'); hold on;%exa10_3 = dsolve('Dy = sin(x*y)', 'y(-1) =0.5', 'x'); %经测试，找不出解析解%ezplot(exa10_3, [-1 1]); %画出解析解的图像%legend('数值解','解析解' )结果：

结果：