ホイン法による常微分方程式の数値解法

を解く。
解析解は、

clear; close all;

x0=0;
xN=1;
N=20;
h=(xN-x0)/N;
x=linspace(x0,xN,N);
y=zeros(1,N);
%syms y0;
a = 1;
y(1)=1;
f = @(x,y) a*y;

for k=2:N
  k1 = h*f(x(k-1),y(k-1));
  k2 = h*f(x(k-1)+h,y(k-1)+k1);
  y(k) = y(k-1) + 1/2*(k1+k2);
end
plot(x,y,'DisplayName','Heun method');
hold on;
plot(x,exp(a*x),'DisplayName','Exact Solution');

legend();

次は、を解いてみる。
解析解は、

clear; close all;

x0=0;
xN=2;
N=10;
h=(xN-x0)/N;
x=linspace(x0,xN,N);
y=zeros(1,N);
a = 1/4;
y(1)=2;
f = @(x,y) -2*a*x*y;

for k=2:N
  k1 = h*f(x(k-1),y(k-1));
  k2 = h*f(x(k-1)+h,y(k-1)+k1);
  y(k) = y(k-1) + 1/2*(k1+k2);
end
plot(x,y,'DisplayName','Heun method');
hold on;
plot(x,2*exp(-x.^2/4),'DisplayName','Exact Solution');

legend();