二阶水箱 matlab 四阶龙格库塔 请问这个二阶常微分方程组用龙格四阶库塔法怎么编写...

要是不用MATLAB自带的ode45函数也可以网上下载一个4阶龙格库塔算法来代替

附上一个仅供参考

function y=DELGKT4_rungekuta(f,h,a,b,y0,varvec) %参数表顺序依次是微分方程组的函数名称，步长，求解范围,初始值向量，变量

format long;

N=(b-a)/h; %求解步数

y=zeros(N+1,1);

y(1)=y0; %赋初值，可以是向量，但是要注意维数

x=a:h:b; %求解范围

var=symvar(f);

for i=2:N+1

K1=Funval(f,varvec,[x(i-1) y(i-1)]);

K2=Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]);

K3=Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K2*h/2]);

K4=Funval(f,varvec,[x(i-1)+h y(i-1)+K3*h]);

y(i)=y(i-1)+h*(K1+2*K2+2*K3+K4)/6;%按照龙格库塔方法进行数值求解

end

format short;

end

***********************

function fv = Funval(f,varvec,varval)

% var=symvar(f);

% if length(var)<4

% if var(1)==varvec(1)

% fv=subs(f,varvec(1),varval(1));

% else

% fv=subs(f,varvec(2),varval(2));

% end

% else

% fv=subs(f,varvec,varval);

var = symvar(f);

varc = symvar(varvec);

s1 = length(var);

s2 = length(varc);

m =floor((s1-1)/3+1);

varv = zeros(1,m);

if s1 ~= s2

for i=0: ((s1-1)/3)

k = strfind(varc,var(3*i+1));

index = (k-1)/3;

varv(i+1) = varval(index+1);

end

fv = subs(f,var,varv);

else

fv = subs(f,varvec,varval);

end

end