Matlab:正则Euler分裂

函数文件1：

function b=F(x0,h,u,tau)
      b(1,1)=x0(1)-u(1);
      b(2,1)=x0(2)-u(2)+2*h*1e8*cos(tau)*x0(2);

函数文件2：

function g=Jacobian(x0,h,tau)
  g(1,1)=1;
  g(1,2)=0;
  g(2,1)=0;
  g(2,2)=1+2*h*1e8*cos(tau); 

函数文件3：

 function x=Euler(h,x0,u,tau)
 % x0 表示迭代初值
 % u 表示上一节点数值
tol=1e-8;
x1=x0-Jacobian(x0,h,tau)\F(x0,h,u,tau);
while (norm(x1-x0,inf)>tol)
    %数值解的2范数是否在误差范围内
    x0=x1;
    x1=x0-Jacobian(x0,h,tau)\F(x0,h,u,tau);
end
x=x1;%不动点

脚本文件：

clear
clc
N=[1e2 1e3 1e4 1e5];
for k=1:length(N)
h=(pi/2)/N(k);
t=0:h:pi/2;
y(1:2,1)=[1;0];%原问题初值
u(1:2,1)=y(1:2,1);%子问题一的初值
for i=1:N(k)
    u(1:2,i+1)=u(1:2,i)+h*[-u(1,i)^2*u(2,i)-u(2,i)^3;u(1,i)+1e8*sin(2*t(i+1))];
    y(1:2,i+1)=Euler(h,u(1:2,i),u(1:2,i+1),t(i+1));
    u(1:2,i+1)=y(1:2,i+1);
end
for i=1:N(k)+1
    Accurate(1:2,i)=[cos(t(i));sin(t(i))];
end
error=abs(Accurate-y);
E_max(k)=max(error(1,:));
end
% plot(t,Accurate(1,:),'-r*')
% hold on
% plot(t,y(1,:),'bp')
% grid on