求解刚性微分方程用什么数值解方法？答案就在这里。|何必手算？Matlab求解微分

clc

close all

clear

syms y1(t) y2(t)

equ1 = diff(y1) == -0.01*y1-99.99*y2;

equ2 = diff(y2) == -100*y2;

equ = [equ1, equ2];

cond = [y1(0)==2, y2(0)==1];

S = dsolve(equ, cond)

fun = @(t, y)[-0.01*y(1)-99.99*y(2); -100*y(2)];

tspan = [0 400];

TimerValue = tic;

[t45, y45] = ode45(fun, tspan, [2, 1]);

timer45 = toc(TimerValue);

plot(t45, y45, 'LineWidth', 2);

hold on;

TimerValue = tic;

[t23tb, y23tb] = ode23t(fun, tspan, [2, 1]);

timer23tb = toc(TimerValue);

plot(t23tb, y23tb, 'LineWidth', 2);

TimerValue = tic;

[t15s, y15s] = ode15s(fun, tspan, [2, 1]);

timer15s = toc(TimerValue);

plot(t15s, y15s, 'LineWidth', 2);

legend('ode45', 'ode23tb', 'ode15s');

errorfun = @(t, y)sqrt(sum((exp(-100*t)-y(2)).^2)/length(t))+sqrt(sum(((exp(-100*t)+exp(-0.01*t))-y(1)).^2)/length(t));

error45 = [timer45 errorfun(t45, y45)]

error23tb = [timer23tb errorfun(t23tb, y23tb)]

error15s = [timer15s errorfun(t15s, y15s)]