# Solving a 6th order non-linear differential equation in Matlab

54 views (last 30 days)
Wissem-Eddine KHATLA on 6 Apr 2022
Hello everyone,
I am trying to solve a high-order non linear differential equation presented right below with : As boundary conditions, I have the following ones : I am trying to apply a shooting method algorithm in order to solve this equation on , I have converted it in a system of 1-st order differential equations as following : with : According to the boundary conditions we thus have : and we have to guess the remaining ones : I have attached to my post, the algorithms that I wrote which are inspired by biblographical researchs : the main code is called "shoot4nl.m" and I have tried to guess some initial values in order to run the algorithm. Also, I have choses a values of very "high" in order to match the conditions...
By doing so, I obtain the following error :
Warning: Matrix is singular to working precision.
> In newtonRaphson2 (line 23)
In shoot4nl (line 14)
Error using newtonRaphson2
Too many iterations
Error in shoot4nl (line 14)
u = newtonRaphson2(@residual,u);
I assume that means that there is a division by 0 that occurs in the process but I remain without any idea to solve this problem...
Could someone help me or bring his mathematical expertise in order to show me some hints ?
Best regards.
##### 2 CommentsShowHide 1 older comment
Wissem-Eddine KHATLA on 7 Apr 2022
@Bruno Luong Thanks for your reply. @Torsten provided a correction to my first attempt by using various in-built functions but it seems like I got an error stating that the "fsolve" function is not correct...

Torsten on 6 Apr 2022
Edited: Torsten on 6 Apr 2022
I think that nobody in this forum will dive into selfmade software if there are well-tested MATLAB alternatives.
So these few lines of code can get you started - backed with professional software:
bc0 = [1 1 1 1];
bc = fsolve(@fun,bc0);
F0 = [1e-5;0;bc(1);bc(2);bc(3);bc(4)]
etaspan = [-10,10];
[eta,F] = ode45(@fun_ode,etaspan,F0);
plot(eta,F(:,1))
function res = fun(bc)
F0 = [1e-5;0;bc(1);bc(2);bc(3);bc(4)];
etaspan = [-10, 0, 10];
[eta,F] = ode45(@fun_ode,etaspan,F0);
res(1) = F(2,2);
res(2) = F(2,4);
res(3) = F(3,1);
res(4) = F(3,2);
end
function dFdeta = fun_ode(eta,F)
dFdeta = [F(2);F(3);F(4);F(5);F(6);(1/9*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3];
end
Wissem-Eddine KHATLA on 10 Apr 2022
Edited: Wissem-Eddine KHATLA on 10 Apr 2022
@Torsten I tried the bvp4c solver using the shooting method result : it doesn't seem to work well with the following code :
% Solving method using bvp4c solver
% Boundary conditions : F'(0) = F'''(0)= F'''''(0) = 0 and F(eta0) = F'(eta0) = F''(eta0) = 0
eta0 = 20;
eta = linspace(0,eta0,100);
yinit = [3.3517;0;0.5079;0;0.1194;0]; % solution of the OCTAVE ode45 integration
solinit = bvpinit(eta,yinit);
sol = bvp4c(@fun_ode, @bcfun, solinit);
eta = sol.x;
F = sol.y(1,:);
figure
plot(eta,F);
function dF = fun_ode(eta,F)
dF=[F(2);
F(3);
F(4);
F(5);
F(6);
(1/9*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3];
end
function res = bcfun(ya,yb)
res =[ya(2); ya(4); ya(6); yb(1); yb(2); yb(3)];
end
Another user suggests me that the problem is not stiff : I am confused since a similar equation is treated in the litterature so I think this is still possible for me to solve it.
Thank you,

Bruno Luong on 9 Apr 2022
Using bvp4c
eta0 = 10;
eta = linspace(0,eta0,100);
yinit = [1e4;0;0;0;0;0];
solinit = bvpinit(eta,yinit);
sol = bvp4c(@fun_ode, @bcfun, solinit);
eta = sol.x;
F = sol.y(1,:);
plot(eta,F); function dF = fun_ode(eta,F)
dF=[F(2);
F(3);
F(4);
F(5);
F(6);
(1/9*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3];
end
function res = bcfun(ya,yb)
res =[ya(2); % F'(0)
ya(4); % F'''(0)
ya(6); % F'''''(0)
yb(1:3)]; % F(eta0), F'(eta0), F''(eta0),
end
Wissem-Eddine KHATLA on 20 Apr 2022
@Torsten F''(-eta0) doesnt equal 1.35. For instance, here is my plot : Compared to the one of the article : 