Solving First-Order Nonlinear Equations Using Interpolation Method in MATLAB

First, the first-order nonlinear differential equation is expanded into an approximate expression concerning the variable Y using the Taylor series expansion method. Then, the number of interpolation points and the number of iterations are determined, the iterative calculation format is written, and the calculation results are displayed.

%% Solution of the first-order nonlinear differential equation using interpolation method %%
 y'(x)=2*y(x)+3*x/(y+10) 
 %% Taylor series expansion regarding y %%
 Y(n)'-2*Y(n)-3*x/(Y(n-1)+10)^2*Y(n)=3*x*(Y(n-1))/(Y(n-1)+10)^2

%%%%%%%%%%%clear all;close all;clc;%%%%%%%%%%%%%%%
 %% Parameter settings
 a1=0;  %% Left boundary point
 b1=1;  %% Right boundary point
 numPoints=60; %% Number of interpolation points
 err1=1.0e-7;  %% Calculation error tolerance
 x1=-cos(pi*(0:1:numPoints)/numPoints); %% Chebyshev interpolation points
 x2=(b1+a1)/2+(b1-a1)/2*x1'; %% Actual interpolation points
 dX0=weifenjuzhen1(x2,3);  %% Derivative matrix
 dX1=dX0(:,:,1); 
 len=length(x2);
 I1=eye(len,len); %% Identity matrix
 y1=x2+1; 
 Max1=999;
 for Loop=1:1:Max1
     y0=y1;
     AA=dX1-2*I1-3*x2./(y0+10).^2;
     B1=3*x2./(y0+10)+3*x2.*y0./(y0+10).^2;
     AA(1,:)=I1(1,:);
     B1(1) = 1;
     y1=AA\B1;
 end
 figure(1);plot(x2,y1);
 disp('Calculation results');

%%%%%%%%%%%%%%%%%%%%function  dY1=weifenjuzhen1(x1,order1)
 %% order1: Derivative order
 n1=length(x1);
dY1=zeros(n1,n1,order1);
b1=weight1(x1); %% Weight function
 for k1=1:1:n1
     for k2=1:1:n1
         if k1 ~= k2
             dY1(k1,k2,1)=(b1(k2)/b1(k1))/(x1(k1)-x1(k2));
         end
     end
 end
 for k1=1:1:n1
     dY1(k1,k1,1)=-sum(dY1(k1,:,1));
 end
 for k1=2:order1
     for k2=1:1:n1
         for k3=1:1:n1
             if k2 ~= k3
                 dY1(k2,k3,k1)=k1*(dY1(k2,k3,k1-1)*dY1(k2,k3))-...
                     dY1(k2,k3,k1-1)/(x1(k2)-x1(k3));
             end
         end
     end
     for k2=1:1:n1
         dY1(k2,k3,k1)=-sum(dY1(k2,:,k1));
     end
 end
end

function  w1=weight1(x1) %% Weight function
 m1=length(x1);
 for k1=1:1:m1
     w1(k1)=1;
     for k2=1:1:m1
         if k1 ~= k2
             w1(k1)=w1(k1)/(x1(k1)-x1(k2));
         end
     end
 end
 w1=w1/abs(w1(1));
end